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PREFACE 


The fundamental concepts and development of computational schemes for the 
solution of parabolic, elliptic, and hyperbolic equations are established in Volume I 
and are extended to the Navier-Stokes equations in Volume II. The primary goal in 
the third volume is to review the fundamentals of turbulence and turbulent flows and 
to extend the governing equations and numerical schemes developed in Volume II 
to include turbulence. 

This volume begins with the basic definitions and concepts in turbulence and 
turbulent flows. Subsequently, the modification of the governing equations and 
numerical schemes is introduced. There are three approaches by which turbulent 
flowfields may be computed. The first approach is based on the averaged Navier- 
Stokes equations either in the form of Reynolds-Averaged Navier-Stokes (RANS) 
equations or Favre-Averaged Navier-Stokes (FANS) equations. These formulations, 
- along with several turbulence models and numerical considerations for the solution 
of equations, are presented in Chapter 21. The second and third approaches are the 
Large Eddy Simulations (LES) and the Direct Numerical Simulations (DNS), which 
are presented in Chapter 23. Since typical computations involved in turbulence and, 
in particular, in DNS, require higher-order schemes such as compact finite difference 
formulations, these formulations are introduced in Chapter 22. 

Finally, a computer code based on the RANS equations and several turbulence 
models have been developed and included in the text Student Guide for CFD- 
Volume III. 

Again, our sincere thanks and appreciation are extended to all individuals ac- 
knowledged in the preface of the first volume. Thank you all very much for your 
friendship and encouragement. 

Klaus A. Hoffmann 
Steve T. Chiang 





Chapter 20 


Turbulence and Turbulent Flows 





20.1 Introductory Remarks 


The vast majority of applications in fluids involves turbulence. Everyday experi- 
ences such as fluid flow in a pipe, flow over an airfoil, flow processes in combustion, 
paint spraying, etc. will exhibit a disorderly complex motion defined as turbu- 
lent flow. Since the majority of fluids involve turbulence, turbulence flow control 
becomes an important aspect of any design process involving fluids. In some appli- 
cations it is advantageous to promote turbulence in order to achieve certain design 
goals. Introduction of dimples on a golf ball is an example whereby enhancing tran- 
sition to turbulent flow reduces flow separation region, and therefore total drag is 
decreased. 

Due to the complexity of turbulent flow and difficulties in its understanding 
and physical interpretation, introduction of a unique conceptual model is difficult 
at best. However, there are certain characteristics and properties associated with 
turbulence which have been shown experimentally and most recently numerically, 
which are used commonly in describing turbulence and turbulent flows. A brief 
description of the physics of turbulence, its generation and development, and a 
conceptual model is presented in this chapter. However, keep in mind that this de- 
scription is not unique and others are provided in literature. The conceptual model 
described in this chapter is provided to shed some light into the world of turbu- 
lence. A review article by Robinson [20.1] is an excellent source for description of 
several conceptual models for turbulence and is highly recommended. Following the 
description of Ref. [20.1], a typical conceptual model is described to identify some 
of the physical characteristics of turbulent flows. Furthermore, several categories 
of turbulent flows are described and some well-established properties of turbulent 
flows are summarized. Finally, the three categories of computational approaches 
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associated with turbulent flows are described. 

A turbulent flow may be defined as a flow which contains self-sustaining fluctu- 
ations of flow properties imposed on the main flow. There are several factors which 
may cause an originally laminar flow to transition to turbulence. The fundamental 
quantity in describing transition to turbulent flow is the Reynolds number. For 
example, for internal pipe flow a transition to turbulent flow can be achieved at 
around a Reynolds number of 2300 and for a boundary layer flow over a flat plate 
at a Reynolds number of around 300,000 ~ 500,000. Obviously, several factors affect 
transition to turbulent flow including freestream turbulence, pressure gradient, heat 
transfer (cooling or heating), surface roughness, and surface curvature. Turbulent 
flows can be broadly categorized into three groups: boundary layers, shear layers, 
and grid-generated turbulent flows. The first two categories are most common and 
will be the focus of our study. 

The first category of turbulent flows is the wall-bounded flows. In these types of 
turbulent flows, the majority of turbulent kinetic energy is produced near the wall 
region. Furthermore, the small scale eddies which are dominant in the near wall 
region are more organized and have similar structure. The wall-bounded turbulent 
flows can be further subcategorized as turbulent boundary layer or fully developed 
turbulent flow. A turbulent boundary layer is simply bounded by a wall and a free 
stream, whereas a fully developed flow is bounded by surfaces, for example, flow in 
a channel or flow in a pipe. 

The second category of turbulent flows is the shear layer. This type of turbulent 
flow grows in the streamwise direction and typically develops universal characteris- 
tics, which may be considered as self-preserving. The turbulent shear layers may be 
subcategorized into three groups as free shear layers, jets, or wakes. A typical mean 
flow velocity profiles for the three categories of shear flows are shown in Figure 20-1. 

The flow field for the flows described above would of course be considerably 
more complex if the flow is supersonic. Such a flowfield will include expansion 
waves, compression waves, shock waves, their interactions with each other, and 
with the shear layer. A typical flowfield for the supersonic wake flow is shown in 
Figure 20-2. 

The third category of turbulent flow is the grid generated turbulent flow, which 
can be produced by passing an initially uniform and irrotational flow through a 
grid composed of rods. The vortices generated by the rods interact together and 
degenerate into turbulence. The grid-generated turbulence is typically isotropic, 
that is, it has a preferred direction. In general, most turbulent flows are anisotropic. 
However, the assumption of isotropy is used in many applications. 
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a. Jet flow b. Free shear layer flow 





c. Wake flow 


Figure 20-1. Illustration of three categories of shear flows. 





The simple turbulent boundary layer over a flat plate and a free shear layer 
will be used to introduce conceptual models. Several processes associated with 
the production of turbulence will be described in this chapter. Before attempting 
to introduce a conceptional model, a brief description of transition is warranted, 
which is given next. Subsequently, some fundamental concepts and definitions are 
reviewed leading to introduction of a conceptual model. 
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Figure 20-2. Schematic of a wake flow for supersonic flow. 





20.2 Fundamental Concepts and Definitions 


Before attempting to describe a conceptual model and mathematical models, 
it is necessary to review several fundamental concepts, physical observations, defi- 
nitions, and terminology related to turbulence and turbulent flows. This goal will 
be accomplished primarily by considering turbulent boundary layers. 

A turbulent boundary layer is commonly divided into several regions, where 
within each region specific turbulence behavior can be identified. A very thin region 
near the surface is referred to as viscous sublayer. The outer region where the flow 
is turbulent is called the fully turbulent zone. A region, which connects these two 
zones, is called the buffer zone. 

The viscous sublayer has been referred to as the laminar sublayer as well. How- 
ever, recent experimental and numerical investigations have shown that in fact this 
region is not laminar. Similarly, the description of buffer zone, which in the past 
has been referred to as a region of transition from laminar to turbulent, is not ac- 
curate. The three regions across a turbulent boundary layer defined above is used 
extensively in description of turbulent flows as well as in development of turbu- 
lence models. Further description of each region based on certain characteristic of 
turbulence development will be presented throughout this section. 

Another categorization of a turbulent boundary layer is to divide it into inner 
and outer regions. The inner region includes the viscous sublayer, buffer zone, and 
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part of the fully turbulent zone. The remaining part of the turbulent boundary 
layer is considered to be in the outer region. 

The division of a turbulent boundary layer to several regions defined above is 
commonly identified by definition of a non-dimensional velocity ut and normal to 
the surface spatial coordinate y*. These quantities are defined as 


u 
ut = — 

Ur 

and i 
+ T 

y —y-— 

v 


where u, is know as the friction velocity given by 


Tw 
Up = |= 
p 


and 7, is the wall shear stress. Now, the various regions are identified as 


yt < 2-8 viscous sublayer (20-1) 
2~8<yt < 2~50= buffer zone (20-2) 
y* > ~50= fully turbulent zone (20-3) 
and 
y* < 100 ~ 400 = inner region (20-4) 
yt > 100 ~ 400 = outer region (20-5) 


Observe that the division between each region in terms of y* is approximate, and 
different values have been used by investigators to identify each region. However, the 
range given above is the most common. The various regions of turbulent boundary 
layer are summarized in Figure 20-3. 


20.3 Introduction to Transition 


Majority of flows begin as orderly fluid motion known as laminar flow. As the 
flow Reynolds number is increased, instabilities within the boundary layer are gen- 
erated. Subsequently, several physical phenomena are developed which eventually 
will lead to transition from laminar flow to a disorderly and random fluid motion 
known as turbulence. Several factors may enhance or delay the transition process. 
Examples of such enhancement mechanism are surface roughness, heat transfer, 
pressure gradient, and freestream turbulence. 
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Figure 20-3. Typical nondimensional velocity profile for a turbulent 
flow over a flat plate. 








In general then, transition is defined as a change where instabilities within a 
laminar boundary layer are initiated. Subsequently several processes such as vor- 
tex stretching, vortex breakdown, and formation of turbulent spots take place, and 
the flow becomes fully turbulent. The entire changes and processes which occur 
over time and space is called transition. Transition process is a complex set of 
events which is extremely difficult to analyze or develop theoretical models for its 
prediction. With the exception of direct numerical simulation (DNS) which is cur- 
rently very limited in applications, no theoretical model for prediction of turbulence 
and its development exists. The only available approach is semi-empirical models 
developed based on experimental data. Our objective at this point is to briefly 
review some of the physical phenomena, which occur during transition. This goal 
is accomplished first by visiting the small disturbance stability theory, because the 
theory can predict onset of instabilities for idealized cases. Subsequently a brief 
description of the processes, which completes the process of transition, is provided. 
We will proceed with description of several relevant stages without detailed mathe- 
matical work. À description of stability theory is given first, followed by description 
of processes, which completes transition. 
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20.3.1 Stability Theory 


In this section the stability theory is reviewed. The mathematical details are 
deleted; instead the procedures to obtain the governing equations and the conclu- 
sions drawn from the solution of the equations are emphasized. The mathematical 
details may be found in texts by White [20.2], Panton [20.3], and Landahl and 
Mello-Christensen [20.4]. 


The governing equation of stability is primarily derived for an incompressible 
flow with constant properties. Assume a vector Q to represent the flow variables and 
Q' to represent perturbations (disturbances) imposed on the flow. The flow variable 
Q + Q' is substituted into the Navier-Stokes equations, and the original equations 
written in terms of Q are subtracted from it. The equations thus obtained involve 
the perturbation variable Q' as the unknown. The resulting equations will be non- 
linear in perturbation properties Q’, thus, typically the equations are linearized by 
neglecting the higher order terms in Q'. This is in fact a reasonable assumption, 
because the disturbances are assumed to be small. Now, if for a given problem 
the disturbance Q' will decay, the problem is stable. On the other hand, if the 
disturbance Q' grows with time, the problem is unstable. The governing equation 
of stability may be solved numerically for a wide range of applications. Tradition- 
ally, analytical solutions of stability equation have been achieved for special cases 
where additional assumptions had to be imposed to reduce the equation. One such 
solution, which is relevant to boundary layer flows, is reviewed at this point. 


The problem of interest is the nearly parallel viscous flow. With this assumption, 
the stability equation is reduced to a single ordinary differential equation. Typically 
a disturbance in the form of an exponential function is introduced and the resulting 
equation is known as the Orr-Sommerfeld equation. 


At this point, the fundamentals and procedures for obtaining the governing 
equations of stability have been reviewed. Again deleting the mathematical details 
or the solution of Orr-Sommerfeld equation for parallel flow, several observations 
regarding the solution of the equations are presented. 


Typical solution of the Orr-Sommerfeld equation is given graphically in Figure 
20-4 for a Blasius boundary layer flow. The horizontal axis represents Res. where 
Res = (puó*)/u, and the vertical axis is a5*, where o is the wave number in the 
z-direction and 6* is the displacement thickness. A curve defining the boundaries of 
the unstable waves which separates the stable and unstable regions is called marginal 
(or neutral) stability curve. A point on the marginal stability curve corresponding 
to lowest Res. is called the critical point. Within the curve of marginal stability is 
the range of unstable wave numbers. These unstable waves are known as Tollmein- 
Schlichting waves, which represent the first indication of instability. 
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Figure 20-4. Graphical presentation of the neutral stability curve for a 
Blasius boundary layer. 





An interesting result obtained from the stability equation applied to Blasius 
boundary layer reported by Squire [20.5] is that when considering a three-dimensional 
disturbance, the corresponding two-dimensional waves are more unstable. This con- 
clusion is important in that two-dimensional analysis can be used to determine the 
limit of stability, which is also valid for the three-dimensional case. 


The following important statement concludes the discussion on stability. The 
stability theory can predict lowest Reynolds number at which instability can be 
initiated. However, note that this is not the location of transition to turbulence. 
Occurrence of instability is only the first indication of a process to transition. In 
fact, typically transition occurs at a distance of 10 to 20 times the distance of 
Zeritical: Furthermore, it should be noted again that several factors affect transition. 
For example, surface roughness, pressure gradient, freestream turbulence, surface 
curvature, and heat transfer have a dominant effect on the process of transition. 
The stability analysis detailed here and a description of transition to follow is for 
a smooth plate with no external factors imposed. The purpose at this point is to 
present a descriptive model for transition in its simplest form. There are still a 
very large number of physical phenomena with regard to transition and turbulence, 
which are not completely understood. 
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20.3.2 Transition 


Based on the discussion of stability just completed, it is well known that at some 
critical Reynolds number instabilities in the form of Tollmein-Schlichting waves 
which are two-dimensional are developed. Shortly thereafter, the unstable waves 
become three dimensional and nonlinear effects come into play. As mentioned pre- 
viously, these instabilities by themselves do not define transition. In fact several 
mechanisms and processes subsequent to initial instabilities must take place to 
complete the process of transition from laminar to turbulent flow. Once the three- 
dimensional instabilities are set in, the nonlinearity effects and interaction with the 
mean flow, result in changes in the mean flow velocity and spanwise stretching of 
vorticity. In addition internal shear layers are formed which are highly unstable. 
Subsequently a breakdown process sets in where the stretched vortices breakdown 
to smaller vortices with random frequencies and amplitudes, producing localized re- 
gions of high shear and intense fluctuations in the shear layer. As a result, random 
and disorderly motion termed turbulence is generated which travels downstream 
and grows in coverage. This process occurs in localized regions surrounded by lam- 
inar flow and is called turbulent spot. Eventually the number of spots is increased 
in the streamwise direction and the flow becomes fully turbulent. This in turn com- 
pletes à process, which is termed transition. This brief description of transition is 
based on a flow over a smooth flat plate. The transition process can be enhanced or 
delayed by several factors identified earlier. A graphical presentation of transition 
process as proposed by White [20.2] is shown in Figure 20-5. 


20.4 Conceptual Model for Turbulent Flows 


A model developed to describe the physics of a problem is called conceptual 
model. A conceptual model for turbulent flow may be developed based on exper- 
imental observations and complemented by numerical observations. At this point 
in time, there are several aspects of turbulent flows, which are poorly understood. 
As a result there is clearly a need for improvement of proposed conceptual mod- 
els. Furthermore, due to uncertainties and the lack of clear understanding and 
interpretation of available data, there are different and sometimes contradictory 
explanations of events related to turbulence. With these remarks in mind, it is 
then obvious that several conceptual models for turbulent flows can be and have 
been proposed. These models are updated as we gain more understanding about 
turbulent flows and it is anticipated that they will eventually converge to a widely 
acceptable model. For the time being though, there are several conceptual models 
proposed by investigators. The objective here is not to review all the conceptual 
models. Instead the goal is to review and familiarize the reader with some of the 
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termologies used in description of turbulence. 

The flow very near the surface in the sublayer and buffer regions moves as streaks 
of low and high-speed fluid. Why these streaks form is still not well understood, 
although Blackwelder [20.6] has suggested that 0.25 jm (1075 in) size dimples (local 
depressions) are sufficient to produce the critical conditions necessary for vortex 
formation. Associated with the low and high-speed streaks are counter rotating 
streamwise vortices, as shown in Figure 20-6. 
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Figure 20-5. Illustration of idealized transition process on a flat plate. 
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Outer region: Large 
scale eddy 


Near wall: Streamwise vortices 
and low speed streaks 


(a) Structure of turbulent boundary layer [20.7] 


Streamwise 
Low speed Vortices 
Streak 





(b) Near wall region has low and high speed streaks and the counter-rotating stream-wise 
vortices are associated with these streaks. [20.7] 


Breakdown 
TN 
oy) 
Oscillation d. ) 


(c) low speed streak lift up off the wall leading to breakdown, known as bursting. [20.7] 


Figure 20-6. Conceptual model for Turbulent boundary layer. 
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Intermittently, the low speed fluid is ejected outward accompanied by inrush of high- 
speed fluid. The ejection of the low speed flow, which is accompanied with oscillation 
and breakup, is termed bursting. It is observed that majority of turbulence energy 
is produced in the sublayer/buffer zone regions [20.8]. 

It is well known that the edge of a turbulent boundary layer is irregular with 
peaks and valleys on the scale of boundary layer thickness. Freestream fluid is 
entrained into the turbulent boundary layer at the irregular interface. Furthermore, 
it is observed that the flow in the outer region of the boundary layer under the bulges 
develops rotational spanwise eddies [20.9]. This outer region large-scale eddy is 
illustrated in Figure 20-6. 

There is ample evidence that the turbulence in the inner region is self-sustaining 
and that the outer layer has some effect on the near wall production process [20.10]. 
Obviously, there is an interaction between the inner and outer regions in that there 
is a transfer of mass, momentum, and energy from the outer layer to the inner layer 
and vice versa. 

Vortical structures have also been used extensively in description of conceptual 
models of turbulent boundary layers. A brief and limited discussion of voritcal 
structures is described in this section. The purpose again is to familiarize the reader 
with some of the proposed models. For a comprehensive review Reference [20.1] 
is highly recommended. The proposed conceptual model populates the turbulent 
boundary layer with various vortices, which are the elements in the turbulence 
production and dissipation as well as transport of mass and momentum between 
the inner and outer layers. One of the most common vortical structures used is 
the hairpin vortex (in high Reynolds number region) and the horseshoe vortex (in 
low Reynolds number region) within the inner region. These vortices are inclined 
in the downstream direction at relatively small angles from the wall, in the range 
of 10° ~ 50°. The formation of these vortices can be visualized as lifting of the 
streamwise vortices, which are stretched and deformed to a horseshoe vortex. As 
discussed previously, spanwise vortices or large eddies are developed in the outer 
region. 


20.5 Production, Diffusion, and Dissipation of Turbu- 
lence 


Any proposed mathematical model must be based on correct representation 
of the physical processes involved in the problem under consideration. For exam- 
ple, the Navier-Stokes equations include unsteady, convective, and diffusive terms. 
Similarly in developing a turbulence model, physical processes involved must be 
identified and must be included in the proposed mathematical model. These pro- 
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cesses include unsteadiness, convection, production, diffusion, and dissipation of 
turbulence. A basic description of production, diffusion, and dissipation of turbu- 
lence is to follow. 

Turbulence is produced by different processes depending on the physics of the 
problem. For example, turbulence is produced in the inner region of a turbulent 
boundary layer by a process called bursting or one may argue the production is the 
result of formation of hairpin vortices. The turbulent eddies formed in the inner 
region are typically of small scale. On the other hand, turbulence production in 
the outer region and formation of eddies are of larger scales. Processes also occurs 
where smaller eddies may grow and become larger eddies and vice versa where the 
larger eddies become smaller and smaller. Thus, one aspect of any turbulent flow 
is production of turbulence, which is commonly referred to as production. 

'The second aspect is diffusion. Recall that mass, energy, and momentum are 
transported within a fluid domain due to random motion of molecules, given rise 
to mass diffusion, thermal diffusion, and viscous diffusion. It can be argued that 
due to the random motion of eddies in a turbulent flow, there is a transport of 
fluid properties from one region of the flow to another region. This process tends 
to increase the mixing process in the fluid and is referred to as eddy diffusion or 
simply as diffusion. Finally, recall that eddies of different scales are formed in a 
turbulent flow. The largest scale in a turbulent boundary layer has an upper limit of 
boundary layer thickness. Similarly, there is a lower limit beyond which the small- 
scale eddy is self-destructive. Of course, molecular viscosity is the primary factor. 
Note that as smaller and smaller eddies are formed, due to their large velocity 
gradient, the viscous forces become considerable and eventually the smallest eddies 
are destroyed. Furthermore, due to the large velocity gradient associated with 
smaller eddies; they dominate the energy dissipation process in the flow. This 
process of energy dissipation and destruction of small eddies is known as dissipation. 


20.6 Length and Time Scales of Turbulence 


It is obvious at this point that there is a wide range of scales associated with 
turbulent flows. The small and large eddies developed in a turbulent flow each 
possess certain characteristics which are important if one attempts to develop pre- 
diction methods for turbulent flows. For example, in direct numerical simulations 
it is imperative that the size of spatial grid be sufficiently fine in order to accurately 
capture the essence of small scale eddies. 

Recall that most of the energy in a turbulent flow is contained in large eddies 
due to their intense velocity fluctuations. Furthermore, the lifespan of the large- 
scale fluctuations is relatively long. On the other hand, the energy associated with 
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smaller eddies is smaller due to their smaller velocity fluctuations and their life 
span is considerably shorter. Furthermore, recall that kinetic energy is transferred 
from larger eddies to smaller and smaller eddies the so called cascading process. 
Eventually, the kinetic energy is dissipated into heat due to the action of molecular 
viscosity at the scale of smallest eddies. 


Now consider a turbulent boundary layer with a freestream velocity of Up and 
a boundary layer thickness ô. Define a turnover time for the large eddies by 6/Uo. 
The energy of the large eddy is proportional to Ug. Hence, the rate of dissipation 
€ would be proportional to U3/6. 

The scales of turbulence can be easily established for smallest eddies if one uses 
Kolmogorov universal equilibrium theorem which states that the rate of transfer of 
energy from larger eddies to smaller eddies is approximately equal to the dissipation 
of energy to heat by the smallest eddies. Therefore, the primary parameters for 
smallest eddies are the mean dissipation and kinematics viscosity. Hence, one can 
use dimensional argument to establish the following length (n), velocity (v), and 


time (T) scales, 
| p> 1/4 
"T Xs 


v = (ve) 


y 2 
ic 
€ 


which are known as the Kolmogorov scales of length, velocity, and time. The 
Kolmogorov length scale 7 is then a measure of smallest eddies which can withstand 
the damping effect of molecular viscosity. It is important to note that the smallest 
length scale of turbulence that we have defined as the Kolmogorov length scale is 
still much larger than the mean free path of molecules, the primary requirement of 
continuum model. Therefore, the continuum model is commonly used for turbulent 
flows. 


It is appropriate at this time to review the concept of isotropy, which is used 
extensively in turbulent flows, and it has significant implications for turbulence 
models and large eddy simulation. Isotropy generically implies invariance with 
respect to orientation. The assumption of local isotropy is used for small turbulence 
scales. Experimental and theoretical considerations indicate that this assumption 
is reasonable under certain circumstances. Nonetheless, the assumption of local 
isotropy is used extensively. The implication of this assumption on the accuracy of 
solutions and alternate procedures are currently under intense investigations. 
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Figure 20-7. Development of Kelvin-Helmholtz vortex. 





20.7 Free Shear Layer Flows 


The general characteristics of turbulent free shear layers are reviewed in this 
section. Consider two flows at different velocities at the upper and lower surfaces 
of a splitter plate. Once the flows pass the trailing edge, instabilities known as 
Kelvin-Helmholtz instability is developed at the interface. Skematic of this process 
and formation of Kelvin-Helmholtz vortex is shown in Figure 20-7. Based on the 
stability analysis discussed previously, it can be shown that the primary factor 
in the development of these instabilities is governed by inviscid character of the 
flow. That is, molecular viscosity has little influence in the development of these 
instabilities as long as molecular viscosity is small Recall that on the contrary, 
the primary factor in the development of wall bounded instability was molecular 
viscosity. Once the Kelvin-Helmholtz instability is initiated, due to the roll up and 
pairing, Kelvin-Helmholtz vortices are formed which are carried downstream. These 
large scale eddies are very regular and possess certain degree of spatial organization. 
Due to their regular formation, they are referred to as coherent structures. These 
large coherent structured eddies which encompass the shear layer thickness, will 
grow linearly as they move downstream due to entrainment of external flow and 
pairing process and swallowing of existing vortices in the shear layer. It is observed 
experimentally that subsequent to initial pairing of vortices, a three-dimensional 
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breakdown into turbulence takes place. Thus, imposed on the quasi two-dimensional 
coherent eddies are the three-dimensional small scale eddies. As in the case of 
wall-bounded turbulent flows, a significant fraction of turbulent kinetic energy and 
Reynolds stresses are due to large scale eddies. 


20.8 Numerical Techniques for Turbulent Flows 


The conservation equations in differential form referred to as the Navier-Stokes 
equations can be solved numerically to predict transition and evolution of turbu- 
lence. In fact the Navier-Stokes equations can be numerically solved for any turbu- 
lent flow and it is common to refer to the approach as Direct Numerical Simulation 
(DNS). This sounds too good to be true and in fact it is at this time. Eventually, as 
more powerful computers capable of massive computations and having capacity of 
huge storage are developed, this dream will become a reality. When this will occur 
is an open question, perhaps within the next 50-4100 years! Let's be optimistic. 

Several issues must be considered and addressed if one is to use DNS for turbu- 
lent flows. First, all scales of turbulence from smallest to largest must be accom- 
modated. This is where limitations to DNS are imposed by present day computers. 
To identify the problem, recall that the smallest scale of turbulence to consider is 
the Kolmogorov length scale 7, and one can consider the boundary layer thickness 
6 as the largest scale. To adequately resolve small-scale eddies whose size is of the 
order of Kolmogorov length scale, a minimum of 4 to 6 grid points are required in 
each direction. It is rather simple to establish an estimate of the number of grid 
points required for a DNS of turbulent flows. For example, consider a turbulent 
boundary layer where the largest length scale is the boundary layer thickness 6. 
For a three-dimensional computation, the number of grid points with a domain of 
6 x 6 x 6 would be approximately 


(rm 


where 6 grid points in each direction are used to resolve small scale eddies. Substi- 
tution of n = (v3/e)'/4 and e = U$/6 yields 


(2) Eco" EG] - 9" une 


Res = D. 
v 





where 


Thus N is proportional to (Res)? for the three-dimensional DNS of turbulent 
boundary layer. Simple estimate of the number of grid points for Res = 5000 
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would be around 4.5 x 10'° grid points. Observe that this is the number of grid 
points required for a cubical domain of 6 x 6 x 6. However, even for a simple 
problem of flat plate or a channel flow, the domain of solution must be in the 
order of several 6 in each direction, most likely several hundred 6 in the streamwise 
direction. And of course for complex geometries it would be in the order of several 
million. It is thus seen that the number of required grid points for relatively complex 
geometries at realistic Reynolds numbers are currently beyond the capability of 
available computers. Current applications of DNS have been limited to simple 
geometries at relatively low Reynolds numbers. A review of these applications will 
be presented in Chapter 23. We can then conclude that, due to the above-mentioned 
limitations, application of DNS to a majority of problems of interest is ruled out 
at present. However, this conclusion does not make DNS obsolete, on the contrary, 
the available results from DNS provide valuable data on the process of transition 
and development of turbulent flow. These data can be used for development and 
improvements of turbulence models. Furthermore, DNS provides wide range of 
data, some of which is very difficult or impossible to measure experimentally. Thus, 
DNS plays an important role in our understanding and prediction of turbulence and 
its importance will increase as more powerful computers are developed. 


In addition to computer requirements of DNS, several numerical issues must 
also be addressed. First, the numerical scheme used must be higher order accurate. 
Recall that typical schemes used in CFD applications are second-order accurate 
in space. In order to accurately resolve various scales of turbulence an order of 
accuracy of four or preferably six is required. For this purpose, compact finite 
difference formulations appear to be very attractive. The compact finite difference 
formulations will be discussed in Chapter 22. 


The second numerical issue is development of accurate boundary conditions. 
Treatment of boundary conditions will be addressed in Chapter 15. 


At this point in time application of DNS for practical problems and its use in 
design process is ruled out. However, over the last several decades approximate 
procedures have been developed which allow us to solve turbulent flow fields. The 
scheme is based on averaging of the fluid properties whereby the high frequency 
(small-scale) fluctuations are removed. These fluctuating terms are then related 
to the mean flow properties by relations, which are known as turbulence models. 
Turbulence models are primarily developed based on experimental data obtained 
from relatively simple flows under controlled environment. That in turn limits the 
range of applicability of turbulence models. That is, no single turbulence model 
is capable of providing accurate solution over a wide range of flow conditions and 
geometries. Nonetheless, turbulence models are a practical approach for solution of 
turbulent flows and are used extensively. Discussion of several turbulence models 
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and numerical considerations are provided in Chapter 21. 

Finally, a third approach is the large eddy simulation (LES), where large scale 
structure in the flow are directly simulated whereas small scales are filtered out 
and are computed by turbulence models called subgrid scale models. Recall that 
previously it was pointed out that the small scale eddies are more uniform and 
have more or less common characteristics. Therefore, modeling of small scale tur- 
bulence appears more appropriate and the models should apply over wide range of 
applications. From numerical point of view, since the small-scale turbulence is now 
modeled, the grid spacing could be much larger than Kolmogorov length scale. This 
in turn allows applications of LES to larger Reynolds numbers possible. Further 
discussions of LES and DNS will be provided in Chapter 23. 


20.9 Concluding Remarks 


e Small-scale eddies appear to be more organized and have similar structure 
and characteristics in all turbulent flows. 


e Large-scale eddies in shear layers are more or less regular and possess co- 
herent structure even though it is possible for them to become irregular. 
Imposed upon the large scale coherent quasi two-dimensional structures are 
three-dimensional small scale eddies. 


e The growth of large-scale eddies are primarily governed by inertia and pressure 
effects, whereas viscosity is the primary factor for small scale eddies. In fact, 
viscosity is the dominant factor limiting the size of small scale eddies. 


e Turbulence is an unsteady phenomenon and it is inherently three-dimensional. 
Eddies overlap in space and larger eddies carry smaller eddies. 


e The majority of turbulence production in a boundary layer takes place in the 
near wall buffer zone where bursting process occurs. 


e Kinetic energy contained in large eddies are transferred to smaller eddies and 
ultimately dissipated into heat at the level of smallest eddies due to molecular 
viscosity. 


e The transfer of energy from larger eddies to smaller and smaller eddies is 
defined as cascading. 


e Kolmogorov length scale is the smallest length scale of turbulence. It is estab- 
lished based on the equilibrium of rate of transfer of energy from the larger 
eddies to smaller eddies and the dissipation of energy to heat by smallest 
eddies. The Kolmogrov length scale for a typical flow is about 0.1~1mm. 
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e As the flow Reynolds number is increased, the range of eddy sizes becomes 
wider with the smallest eddies getting smaller. 


* An increase in the flow Reynolds number results in longer energy cascade. 


e It is perfectly reasonable and logical to solve the governing equations of fluid 
motion subject to well defined initial and boundary conditions by numerical 
schemes for turbulent flows. This approach is known as direct numerical 
simulation (DNS). All relevant scales of turbulence must be accommodated 
for DNS of turbulent flows. As a result, DNS is extremely computer intensive 
and for most applications beyond the capabilities of available computers. 


e The number of grid points for DNS of turbulent flows scales as Re?/4 for 
three-dimensional problems and as Re? for two-dimensional problems. When 
the smallest scales of the flow are filtered out and the remaining scales are 
directly computed by the filtered Navier-Stokes equation, the procedure is 
called large eddy simulation (LES). The filtered small scales are determined 
with the subgrid scale models and thus introduce some degree of uncertainties 
into the computations. 





Chapter 21 


| Turbulent Flow and Turbulence Models 


21.1 Introductory Remarks 


Fundamental concepts, basic definitions, conceptual models, and numerical 
considerations for turbulent flows were established in the previous chapter. Three 
different techniques for the computation of turbulent flows were identified, and a 
brief description of each was presented. It was pointed out that, due to large com- 
putational time and computational resource requirements for DNS and, to some 
extent, for LES, these techniques are used more or less for research-oriented appli- 
cations at this time. For routine computations involving turbulence, the averaged 
Navier-Stokes equations complemented with turbulence models are utilized. The 
objectives of this chapter are (a) to review the relevant averaged Navier-Stokes 
equations and turbulence models, (b) to explore the appropriate numerical schemes 
for solutions, and (c) to investigate several applications. Discussion of LES and 
DNS are presented in Chapter 23. 


21.2 Fundamental Concepts 


In this section, the fundamentals of turbulent flows, some of which were in- 
troduced in the previous chapter, are explored and reiterated. This explanation is 
accomplished by considering the classical turbulent boundary layers. 

In many high Reynolds number flows, the effect of viscous forces is dominant 
in the regions adjacent to the surface. This region is known as the boundary layer. 
As it was pointed out in the previous chapter, a boundary layer usually starts 
with well-behaved streamlines where the fluid mixing is at a microscopic level. 
This type of boundary layer flow is identified as a laminar boundary layer. Due 
to conditions imposed by the geometry and flow field, such as surface roughness, 
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surface temperature, surface injection, and pressure gradient, the mixing of the 
fluid increases and takes place at the macroscopic level; streamlines are no longer 
well-behaved. This type of flow is known as turbulent flow. There is a transitional 
region between laminar and turbulent boundary layers known as the transition 
region, which was described in Section 20.3. The various boundary layer regions 
are shown in Figure 21-1. 


S 


L— 34. 


Laminar Transition Turbulent 





Figure 21-1 Various flow regimes near the surface. 


Turbulent 


. Figure 21-2 Comparison of typical velocity profiles in laminar and 
turbulent boundary layers. 
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As a result of heavy mixing of the fluid in a turbulent boundary layer and the 
associated large momentum flux, the velocity profile becomes fuller in a turbulent 
boundary layer as compared to a laminar boundary layer; i.e., the velocity gradient 
near the wall for a turbulent boundary layer is larger than the velocity gradient 
of a laminar boundary layer. Typical velocity profiles for laminar and turbulent 
boundary layers are shown in Figure 21-2. 

In addition, boundary layer flows can be classified into velocity boundary layers 
and thermal boundary layers. For most problems, velocity boundary layer thick- 
ness and thermal boundary layer thickness are not identical. Typical velocity and 
thermal boundary layers are shown in Figure 21-3. 


Thermal Boundary Layer 


Velocity Boundary Layer 





Figure 21-3 A typical comparison of velocity and thermal boundary layers. 


Recall that à nondimensional parameter was previously defined as the Prandtl 


number and expressed as 
Pr => =f% 
perl 
a k 


This parameter represents the ratio of momentum transfer in the flow to that of 
heat transfer. Therefore, it may represent relative velocity and thermal boundary 
layer thicknesses. Figure 21-4 illustrates this relation. If 6 and 6; are used to denote 
the velocity and thermal boundary layer thicknesses, respectively, then, 


for Pr«1 ; 626 
Pr=1 ; é= 6 


and for Pr>1 A é< ô 
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Figure 21-4. The velocity and thermal boundary layer thicknesses for 
various Prandtl numbers. 








21.2.1 Universal Velocity Distribution 


A turbulent boundary layer over a smooth flat plate with zero pressure gra- 
dient is perhaps one of the simplest type of turbulent flows. However, even this 
"simple" turbulent flow is in fact very complex, and an analytical solution can not 
be obtained. It is therefore desirable to develop semi-empirical relations, which 
càn be used to provide approximate solutions for the simple turbulent flows. One 
such relation is the Universal Velocity Distribution, which is reviewed in this sec- 
tion. This relation is developed based on dimensional analysis and incorporation of 
experimental data to determine the constants appearing in the expression. 

Physical parameters which influence the velocity profile near the wall include 
density (p), viscosity (u), distance to the wall (y), surface roughness (k), and the 
velocity gradient at the wall [(du/dy),], or, equivalently, the shear stress at the wall 
(tw). Thus, a functional relation for the near wall velocity is written as 


u= F(p, Wy k, Tw) (21-1) 


Note that the influence of the outer layer on the near wall velocity profile enters 
through 


ut = f(y*, k*) (21-2) 


where 
kt = ut (21-3) 


24 Chapter 21 





Expression (21-2) is reduced to 


ut = f(y") (21-4) 


for a smooth wall. This relation is known as the law of the wall. In terms of the 
different regions of turbulent boundary layer defined previously in Section 20.2, the 
law of the wall is valid for the viscous layer, the buffer zone, and the fully turbulent 
portions of the boundary layer. It should be noted that the total shear stress 
composed of viscous stress and Reynolds stress is approximately constant within 
the regime governed by the law of the wall. 

The outer region velocity profile can be expressed as 


U— Ue = G(p, Y, 6, Ur, dp/dz) (21-5) 


or in terms of nondimensional parameters as 
Un Ue ( up 2) (21-6) 


This relation is known as the defect law. Observe that viscosity does not appear 
in the functional relation (21-5). Recall that as previously discussed in Section 20- 
20, the effect of viscosity is to dissipate the small eddies due to existence of large 
velocity gradient. This effect is dominant near the wall where the Reynolds number 
is relatively low. However, the contribution of small eddies to the Reynolds stress 
in the high Reynolds number outer region is in fact negligible. As a consequence it 
is argued that the outer region of the velocity profile is essentially independent of 
viscosity. 





At this point, a relation for the near wall region known as the law of the wall, 
given by the functional form (21-1) or equivalently by (21-2), and a relation for the 
outer region known as the defect law, given by (21-5) or (21-6), has been established 
based on dimensional analysis and experimental observations. It is then obvious 
that these relations must merge together smoothly over a finite region between 
the near wall region and the outer region. This region is known as the log layer 
or the overlap layer. The length of the log layer is problem dependent, and, in 
particular, a function of Reynolds number (u,6)/v. Typical range of log layer is 
~ 50 < yt <~ 400, or, in terms of y/ó, ~ 0.01 < y/ó <~ 0.2. 

Now consider a smooth flat plate with no imposed pressure gradient. Then, 
from relation (21-4), the law of the wall is 


uly (=) =f (5 zt) (21-7) 


Ur V 
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and, from relation (21-6), the defect law is written as 


— ($) (21-8) 


Ur Ur 
Since, in the log layer relations, (21-7) and (21-8) must merge together, we set 
u Urd Ue y 
Z = z 21- 
( 3 v Ha) CES 


Functional analysis indicates the validity of relation (21-9) if both f and g have 
logarithmic form. Thus an equation in the following form is established 





=< tny +B (21-10) 


ur 


where B and «x (von Karman Constant) are nondimensional constants determined 
from experimental data, with typical values of x = 0.4 ~ 0.41 and B = 5.0 ~ 5.5. 
Thus, Equation (21-10) is written as 


ut = 2.5ényt + 5.5 (21-11) 


The equation for the outer layer can be expressed as 


U-Ue l, Y 
n cw A 5 A (21-12) 
where 
A & 2.35 


It is important to emphasize that both relations given by (21-11) and (21-12) 
are valid only for incompressible flows over smooth flat plates with zero pressure 
gradient. 

Now consider the velocity profile very near the wall within the viscous sublayer, 
that is, yt < 2 ~ 8. This is a region where viscous shear dominates and the velocity 
profile can be approximated to be linear. Therefore, 


NES ST 
pa T, u 
fo Hw ey 
p Py 
or s 
MM 
from which 





26 Chapter 21 





or 
ut =yt (21-13) 


The velocity profile in the sublayer region given by (21-13) must smoothly merge 
with the velocity profile for the log layer given by (21-11). A relation satisfying this 
requirement is established using experimental data as 


ut = 5.0£ny* — 3.05 (21-14) 


which is valid for the buffer zone, that is, 2 ^4 8 < y* < 2 ~ 50. 

Now that several relations which essentially form the law of the wall have been 
established, Figure 20-3 is repeated as Figure 21-5, where the equation for each 
region is identified. The compressibility law of the wall for the log layer can be 
derived for most turbulence models. The main assumption used in the derivation 
is based on the fact that the convection effects are small within the log layer, and 
therefore they are ignored in the equations of fluid motion and the turbulence models 
are considered. The mathematical details are omitted here, however the details can 
be found in References (21.1, 21.2). The resulting expressions from the k — € model 
and the k — w model are, respectively, 


ee eee (21-15) 
Ke 
where 
" 25A 
B, - B4 —fn ( ) (21-16) 
Ke Pw 
and 
ut = > tny* +B (21-17) 
where 
i QNA 
B, — B — tn (2) (21-18) 
Kw w 


and p, is the density at the wall. The constants ke and x, are known as the effective 
von Karman constants. In fact, Ke and Kuy are not constants, but are variables for 
which expressions can be developed, as in Reference (21.1). However, the variation 
in the values of x. and «u is typically small, in the range of about 0.2 percent for 
Ku and 5 percent for x, when the Mach number range is up to about 5. Therefore, 
for simplicity, one may select Ku = Ke = K = 0.41. 


Turbulent Flow and Turbulence Models 27 


i — ——————————— Law of the wall region _ se d 
Log layer 


1 
-= Defect law region — 


paee Inner region —— > Outer region 


| 


|æ Viscous pæ Buffer zone Exe Fully 
sublayer turbulent 
20 + 





Eq. (21-13) 
15 1 
Eq. (21-11) 











Figure 21-5. Nondimensional velocity profile for an incompressible turbu- 
lent flow over a flat plate and identification of different regions 
within the turbulent boundary layer. 





21.8 Modifications of the Equations of Fluid Motion 


In order to include and account for the effect of turbulence in a flow field, the 
equations of fluid motion are modified and amended by turbulence models. There 
are two approaches to reformulate the Navier-Stokes equations for this purpose. In 
both approaches, an averaging process is used. The resulting equations are known 
as the Reynolds Averaged Navier-Stokes equations (RANS) and the Favre Averaged 
Navier-Stokes equations (FANS). The governing equations and assumptions for each 
approach are presented in the section (21.3.1) for RANS and in the section (21.8) 
for FANS. 
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21.3.1 Reynolds Averaged Navier-Stokes Equations 


Traditionally this modification is accomplished by representing the instanta- 
neous flow quantities by the sum of a mean value (denoted by a bar over the variable) 
and a time-dependent fluctuating value (denoted by a prime). Mathematically, it 
is expressed as 


ffe (21-19) 

where "m 
Tod T di 21-20 
t=, f (21-20) 


These quantities are shown graphically in Figure 21-6. The time averaging of à 
fluctuating quantity over a time interval At results in 


tot At 


pot ‘dt ox i 
Pax J f'dt=0 (21-21) 


Figure 21-6. Illustration of the average and fluctuating quantities for 
steady and unsteady flows. 





The time interval, At, used in the definitions above must be larger than the period 
of fluctuating quantities but smaller than the time interval associated with the 
unsteady flow. Thus the time interval is problem dependent, i.e., depends on the 
geometry and physics of the flow field being investigated. Note that for a steady 
flow the time-averaged mean value is constant while for an unsteady flow it is a 
function of time. 

In order to carry out the mathematical details, the following averaging rules are 
applied: 


fj =f (21-22a) 
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f+9 = f*9 (21-22b) 
fg = fs (21-22c) 

i - T (21-22d) 

F #0 (21-226) 

fg # 0 (21-22f) 
G+" = P+? (21-228) 


It is important to recognize that the averaged values of the products of fluctuat- 
ing terms represented by f'g’ are not, in general, zero. Indeed, these quantities—in 
particular, those involving the velocity fluctuations—are critically important and 
represent the cornerstone of turbulent flow effects. 

The conservation of mass for a two-dimensional flow is given by the first com- 
ponent of Equation (11-192), which can be expressed in dimensional form as 

op 


FE + Zio + 5-00) + (po) =0 


The density and velocity components are replaced by relation (21-19) to provide 
fa] 1 0 fx t ü — t / 
a 0 A) + PHP) Go D e i Gg) Gr) 


typ dev) =0 


Once this equation is expanded, it is time averaged according to rules set by 
(21-22) to provide 
db + E (put Fu) + 5 (pe 70) + (pvo 77) =0 (21-23) 
The mathematical manipulation for the momentum equation is considered for 
the z-component, and, subsequently, the results are extended to the y-component. 
The z-component of the momentum equation in dimensional form is given by the 
second component of (11-192) as 


9 EA 24 Eur jee am 40u 28v 
gi Pt gu TP ay rr y ^" = Or |" V 38x 3 Oy 


4:9. | [96 ,.0vY| | LAELA -2 B f t 
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The instantaneous values are replaced by time-averaged mean and fluctuating values 


to provide 
3 ^ fz I Ot ^ fs; n2 = ! 
3 9 6) (ie D] sz [o 6) Gi Y 7] 


^ 5, (p+ p) (+u) +v) 4 dd (P+ p) (a+) @+v)] = 


(een ibl) aero 


a On y PEE 2 ð jtn , 
e 


Since the fluctuating component of viscosity is usually small, it has not been in- 
cluded in this equation. Now the entire equation is time averaged and subsequently 
rearranged as 


2 (put ze) (pu par pu + pu + niga) 

O TIT OK a mUS aq A odu + ou 
a, (PUR + Pu + upv + vplu! + p'u'v id e ^v! + up v! + vp! ! 4 gu) = 
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Writing this equation in terms of shear stresses, one obtains 


a (pu+ za) +e (pu? -- p4- up) + 5, (ouv u7) + 





Ó 
(exp — Bu? — pf — pu) + 5 (ra, - vuv - Fay — Up) + 


3 -— =a 
Oz 

al. 2 8 nia — 7 _ Wa e 
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Similarly, the y-component of the momentum equation becomes 
Q Ó (qo | [7 
aL (P+ 9v) + s; (PUT vov) + dy (pv? tp up) 


+2 (pego) = 2 (ra -E - ga - p) 
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The energy equation given by the fourth component of (11-192) in terms of the 
total energy per unit mass is expressed as 


3, (pe) + lioe D ul + pee p) o] S [oe +p) = 


8 à 
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Note that the total energy is assumed to be composed of internal energy and ki- 
netic energy. The mathematical procedure used to include turbulence is similar to 
previous manipulations. The final form of the energy equation becomes 


ð 1l ;/— 1 —— zt 
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Now it is obvious that the number of unknowns in the system of equations 
composed of (21-23) through (21-26) has increased significantly. Before considering 
additional relations to close the system, some reduction of terms is introduced. It 
has been shown experimentally that, for flows up to about Mach five, turbulence 
structure remains unchanged and is similar to that of incompressible flows (the 
so-called Morkovin hypothesis). Therefore, the fluctuating correlations involving 
density fluctuations are usually neglected. In addition, triple correlations such as 
v® are assumed smaller than the double correlations and are neglected as well. 
Finally, Equations (21-23) through (21-26) are expressed as 
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(21-28d) 
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Additional unknowns appear in the system of equations given by (21-27) as a 
result of modification of the original system (11-192) to include the effect of turbu- 
lence. It is obvious that the solution of Equation (21-27), which includes turbulence 
quantities, would be a very challenging task. In fact, in practice, the effect of tur- 
bulence in most applications is simulated by the introduction of turbulent viscosity 
and turbulent conductivity, as will be shown in Section 21.3.1.1. In either case, 
there are additional unknowns in the Navier-Stokes equations, due to turbulence. 

In order to close the system, supplementary relations must be introduced. These 
relations are known as turbulence models, which in general relate the fluctuating 
correlations to the mean flow properties by means of empirical constants. Turbu- 
lence models vary in degree of sophistication from simple algebraic equations to 
systems of partial differential equations. Surprisingly enough, in many instances 
simple algebraic models provide as good a result as sophisticated models, with less 
computation time. Obviously, the ultimate goal is to simulate and resolve all scales 
of turbulence directly by the Navier-Stokes equations. However, due to limitations 
in the capacity of current computers, the direct simulation approach is more or less 
at the research level, and it has not matured for routine computations as yet. With 
the advancement in computer technology, however, this goal will be accomplished 
in the future. Research toward this goal is underway. 


21.3.1.1 Turbulent Shear Stress and Heat Flux: 

In this section, some of the fundamental concepts introduced by the early inves- 
tigating pioneers are explored. For the sake of simplicity, the steady incompressible 
boundary layer equations are considered. The governing equations are expressed as 
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In the equation above, o is the thermal diffusivity. 

It can be shown that the macroscopic momentum exchange due to turbulence is 
represented by — pu'v', which is referred to as the turbulent shear stress (or Reynolds 
stress), denoted herein by 7. It is customary to combine the laminar and turbulent 
Shear stress terms in Equation (21-31) as 

eu 8 m 18 
vay = ay (wv) = d 5, + Ti) (21-33) 


For a laminar flow under the assumptions stated, one may write 


en ape 
In order to express the turbulent shear in a similar form, consider the concept of 
eddy viscosity, where the turbulent shear stress is related to the gradient of the 


mean flow velocity. This analogy was introduced by Boussinesq and is referred to 
as the Boussinesq assumption. Thus, one writes 
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where p (or v) is known as the turbulent viscosity or eddy viscosity. Hence, 
Equation (21-31) may be expressed as 
Ou " 9" l êp __ 1p Op | 
p ðy $ p 
In a similar fashion, turbulent PLOTS is defined to combine laminar and 


turbulent heat fluxes. For a laminar flow, the Fourier heat conduction law is ex- 
pressed as 
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where k; is the turbulent conductivity and œ is the turbulent diffusivity. Now, the 
energy equation may be rearranged as 
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or in terms of turbulent parameters, 
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Note that, up to this point, the problem of closure has not been addressed, 
ie., two additional unknowns still remain. Allthat has been accomplished is the 
rewriting of the equations by introducing the Boussinesq assumption. 

It is well known that in a laminar flow the mixing of fluid is at the molecular 
level, and the viscous stresses and heat fluxes are due to momentum and energy 
transfer by molecules traveling a distance of mean-free path before collision. A 
similar concept is extended to turbulent flows, where it is assumed that lumps of 
fluid travel a finite distance before collision and before losing their identity. The 
resulting momentum and energy exchange give rise to what was defined as turbulent 
shear stress and turbulent heat flux. This finite distance is defined as the mixing 
length, lm. This concept was introduced by Prandtl and is known as the Prandtl 
mixing length. The Prandtl hypothesis is expressed as 
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In terms of turbulent viscosity and turbulent conductivity, the following may be 
written: 


Ou 
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The concept of turbulent viscosity and turbulent conductivity can be easily ex- 
tended to general flow fields governed by the Navier-Stokes equations. For example, 


for a 2-D problem, 
1 
2 2]3 
ji pl I) + (z) | (21-45) 


What has been accomplished up to this point is the introduction of new param- 
eters, and the central issue of closure is not complete as yet since one still has two 
additional unknowns, now in the form of p and k;, or, equivalently, lm and le. 
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In order to write expressions for mixing lengths, experimental investigations are 
heavily relied upon; and J, and le are modeled for the various flow regimes, namely 
the inner and outer regions. 

Before proceeding further, it is beneificial to elaborate on the concept of mixing 
length Im and the turbulence length scale l, because both are used in the develop- 
ment of turbulence models. Recall that previously, in Chapter 20, the turbulence 
length scale was defined as a characteristic length associated with the eddy size in a 
turbulent flow. The smallest of the length scales was defined as Kolmogorov length 
scale associated with the smallest eddy. It was also established that the size of the 
length scale within a flow varies substantially from the smallest scale (Kolmogorov 
length scale) to the largest scale (the size of boundary layer or shear layer thick- 
nesses). Furthermore, length scales are flow dependent and would be different for 
different flow conditions. 

The concept of mixing length proposed by Prandtl was developed based on 
the analogy of momentum exchange in laminar and turbulent flows. Thus, mixing 
length is a nonphysical property. The mixing length concept was used by pioneering 
investigators in the development of turbulence models expressed as algebraic rela- 
tions. Later on, when one- and two-equation models were developed, the concept 
of turbulence length scale was incorporated in the expressions for turbulence quan- 
titites. Thus, the mixing length and the turbulence length scale are in fact different 
quantities, and in general the mixing length would not be equal to the turbulence 
length scale. However, for a special case where the ratio of production to dissipation 
remains constant, the length scale would be proportional to the mixing length, that 
is, | =constant : lm. 

In the introduction of turbulence models to follow, the mixing length concept will 
be used primarily for the algebraic turbulence models, whereas turbulence length 
scale is used in the one- and two-equation turbulence models. 


21.3.1.2 Flux Vector RANS Formulation: 


A simple and practical approach to include the effect of turbulence in the Navier- 
Stokes equations is simply to replace the molecular viscosity u in the stress terms 
by (u + 14) and the term (/Pr) in the heat conduction term by (u/Pr + u;/ Pr;). 
The parameter Pr, is called the turbulent Prandtl number, and for air it is generally 
taken to be 0.9 for wall bounded flows and 0.5 for shear layers. 

Returning to the nondimensional Navier-Stokes equations in the computational 
space given by Equation (11-60) and repeated here for convenience, one has 
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where the flux vectors are given by (11-62) through (11-67). Now, to include the 


(21-46) 


38 


Chapter 21 


effect of turbulence, Equation (21-46) is used along with the following definitions 
for the viscous stress and heat conduction terms. 
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where the value of 44; is provided by a turbulence model. 


21.4 Turbulence Models 


A turbulence model is a semi-empirical equation relating the fluctuating cor- 
relation to mean flow variables with various constants provided from experimental 
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investigations. When this equation is expressed as an algebraic equation, it is re- 
ferred to as the zero-equation model. When partial differential equations are used, 
they are referred to as one-equation or two-equation models, depending on the num- 
ber of PDEs utilized. Some models employ ordinary differential equations, in which 
case they are referred to as half-equation models. Finally, it is possible to write a 
partial differential equation directly for each of the turbulence correlations such as 
wv’, u®, etc., in which case they compose a system of PDEs known as the Reynolds 
stress equations. 

In the following sections, some commonly used zero-equation, one-equation, and 
two-equation turbulence models, as a sample of the wide range of models, are inves- 
tigated. Several recently published texts, survey papers, and comparative analysis 
papers on turbulence models are excellent references for those interested in various 
aspects of turbulence models, [21-1] through [21-8]. 


21.4.1 Zero-Equation Turbulence Models 


The zero-equation models are equations wherein the turbulent fluctuating cor- 
relations are related to the mean flow field quantities by algebraic relations. The 
underlying assumption in zero-equation models is that the local rate of production 
of turbulence and the rate of dissipation of turbulence are approximately equal.. 
Furthermore, they do not include the convection of turbulence. Obviously, this is 
contrary to the physics of most flow fields, since the past history of the flow must 
be accounted for. Nevertheless, these models are mathematically simple and their 
incorporation into a numerical code can be accomplished with relative ease. 

Generally, most models employ an inner region/outer region formulation to rep- 
resent mixing length (the subscript m will be dropped from lm in the following). A 
commonly used model utilizes an exponential function (van Driest damping func- 
tion) for the inner region, whereas the outer region is proportional to the boundary 
layer thickness. Mathematically they are expressed as 

li = k(1— ev /^)y (21-56) 
and 
lo = Cô (21-57) 
where « is the von Karman constant (~ 0.41), and A* is a parameter which depends 
on the streamwise pressure gradient. For a zero-pressure gradient flow, it has a value 
of 26. The constant C, in Equation (21-57) is usually assigned a value of 0.08 ~ 0.09, 
whereas 6 is the boundary layer thickness. 


Another formulation commonly used for the outer region turbulent viscosity is 
the Cebeci/Smith model expressed as 


y= ales . (21-58) 
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where o is usually assigned a value of 0.0168 for flows where Reg (momentum 
thickness Reynolds number) is greater than 5000 and 6° is the kinetic displacement 
thickness defined as 





= T5 ( = =) dy for an incompressible flow, and (21-59) 
6 = n Pi she dy for a compressible flow (21-60) 
o PeUe 


Recall that the Reynolds number based on momentum thickness is defined as 


Rep = d (21-61) 


where the momentum thickness Ó is 





0-— D i = ( — z) dy for an incompressible flow, and (21-62) 
"S age ( =) 3 
d= J PS 1 u dy for a compressible flow (21-63) 


The algebraic model described above requires information regarding the bound- 
ary layer thickness and flow properties at the boundary layer edge. When the 
Navier-Stokes equations are being solved, it may be a difficult task to determine 
the boundary layer thickness and the required properties at the edge. That is es- 
pecially the case when flow separation exists within the domain. However, when it 
is necessary to determine the extent of the viscous region within the domain, the 
total enthalpy is usually used. 

A turbulence model which is not written in terms of the boundary layer quanti- 


ties was introduced by Baldwin and Lomax [21-9]. The inner region is approximated 
by 


m = plo (21-64) 
where w is the vorticity defined as 
ðv ðu 
=e 21-65 
Ox Oy ( ) 


and | is given by (21-56). The nondimensional space coordinate yt can be written 


as 
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or 
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where the subscript w indicates the wall quantities. The outer region is approxi- 
mated by 


He = APCopFwakeF kied (21-66) 
where a is assigned a value of 0.0168 (as in the Cebeci/Smith model) and 


22] 





Fake = min aso > CwakeYmax (21-67) 


A typical value for Cwake is 1.0. In Equation (21-67), the following definitions are 
employed, 


Gmax = max (Zi) . (21-68) 


where « is the von Karman constant, and l, the mixing length, is determined by the 
van Driest function given by (21-56). The difference between the absolute values 
of the maximum and minimum velocities within the viscous region is denoted by 
AV. For wall bounded flows, the minimum velocity occurs at the surface where the 
velocity is zero, then 

AV = (u? + v?) 


For shear layer flows, AV is defined as the difference between the maximum velocity 
and the velocity at the ymax location, that is, 


AV = (u? + v!) — (w+ v2 


Fie is the intermittency factor defined as 
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Ymax 


and Ymax is the y location where Gmax occurs. 

The typical values for the Klebanoff constant Cki and Cop are respectively 0.3 
and 1.6 for zero or mild pressure gradient. However, to include the influence of 
pressure gradient, the following expressions as suggested in Reference [21-10] may 
be used: 
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and u, is the friction velocity. The velocity gradient in f* is calculated outside the 
viscous region. Once the Klebanoff constant is evaluated, Ccp is determined from 


= $e MR 
2Ckieb(2 — 3Ckieb + Chleb) 


Finally, the turbulent viscosity distribution across the boundary layer is deter- 
mined from 


Cor (21-71) 


pi = MIN (Hi, He.) (21-72) 
that is, 1 is set to 44, from the wall to a location for which Ha exceeds the value of 
14, at which point pu is set to u,. 

To model the eddy diffusivity, œ, or equivalently the turbulent conductivity, 
the Reynolds analogy may be used. Recall that the Reynolds analogy assumes a 
similarity between the momentum transfer and heat transfer. Therefore, a turbulent 
Prandtl number is defined as 


Pr = eod (21-73) 


For most flows, it is assumed that the turbulent Prandtl number remains con- 
stant across the boundary layer. For air, Pr; = 0.9 for wall bounded flows and 0.5 
for shear layers. Thus, the turbulent conductivity is determined as 

tcp 
ki = = 21-74 
i Pr, ( ) 
where p is provided by the turbulence models just described. 

To initiate the computations, a transition location zy is specified, and the initial 
value of the turbulent viscosity is set equal to zero everywhere within the domain. 
Subsequently, in regions where z > zr, the turbulent viscosity is computed from 
(21-72) and updated after each time step (or iteration). 


21.4.2 One-Equation Turbulence Models 


As seen in the previous section, zero-equation models employ an algebraic 
relation for the eddy viscosity. These models specify length and velocity scales in 
terms of the mean flow, thus implying an equilibrium between mean motion and 
turbulence. 

One-equation models employ a partial differential equation for velocity scale, 
whereas the length scale is specified algebraically. The velocity scale is typically 
written in terms of turbulent kinetic energy k defined as 


ko +07 uf) (21-75) 
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The turbulent viscosity u is now written as 
ji = pki (21-76) 


In the following section, two recently developed one-equation turbulence models 
are reviewed. 


21.4.2.1 Baldwin-Barth One-Equation Turbulence Model: 

The Baldwin-Barth one-equation model is obtained from the standard k-e two- 
equation model. Note that the k-e model has not been introduced as yet, and it will 
be introduced later on in Section 21.4.3.1. However, since the mathematical details 
of the Baldwin-Barth one-equation model will not be considered here, knowledge of 
the k-e model at this point is not necessary. Mathematical details can be found in 
a report by Baldwin and Barth [21-11]. 

The procedure to obtain the Baldwin-Barth one-equation model is to combine 
the k-e two-equation model to provide a single equation in terms of the Turbulence 
Reynolds number Rer, where Rer = k?/ve. The variables in Rer which are k, 
e, and v are, respectively, turbulence kinetic energy given by (21-75), dissipation 
of turbulence given by (21-152), and kinematic viscosity. To account for the near 
wall regions, the turbulence Reynolds number is split into two parts as Rer = 
Rer f (Rer), where f is a damping function such that Rer = Rer for large Rer. 

The transport equation for turbulence Reynolds number Rer is written as 


2 Rer)+ V.V(vRer) = 


(v+ A) V'(vRe;) — (Vw) - V(vRej) + (cafa - co) V vRepP. (21-77) 


The turbulence production based on Boussinesq assumption is provided from 
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The kinematic eddy viscosity v is determined from 
y= c,(vRer)DiD; (21-79) 


where D, and D; are damping functions which extend the validity of the model to 
near wall regions and are given by 
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Various functions and constants appearing in Equation (21-77) are as follows. 
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Equation (21-77) can be written as 
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as obtained from (21-78) for two-dimensional applications. Due to numerical con- 


siderations, Equation (21-83) can be rewritten by manipulation of the diffusion term 
and with the assumption of constant v as follows: 


a Re V-V(vRer) = 2 (v + A) V*(vRe;) — LV - AV (v Rer) 


+ (Cafa — ca) / e, DiD3SvRer (21-85) 


Note that the unknown in either Equation (21-83) or (21-85) is Rer. Once Rer is 
determined, relation (21-79) is used to determine 14. 


21.4.2.1.1 Nondimensional form: Since the governing equations to be solved are 
typically in nondimensional form, the equations for turbulence models are also 
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nondimensionalized. The same nondimensional variables, as defined in Section 11.3, 
are used, namely, 
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Now, the nondimensional form of the Baldwin-Barth one-equation turbulence 


model corresponding to (21-83) or (21-85) is given by the following equations where 
the asterisk has been dropped. 


a Pe V NUS) = g (v 5) hen - g Lv) voe) 
+ (cafa — ca) / e DiDaSvRer (21-86) 
or l 
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21.4.2.1.2 Baldwin-Barth turbulence model in computational space: In order 
to numerically solve either Equation (21-86) or (21-87) for a general domain, a co- 
ordinate transformation is performed. This transformation is similar to the trans- 
formation of the Navier-Stokes equations given in Chapter 11. The expressions for 
the Cartesian derivatives are 
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Now, for simplicity, define 


vRer = R 


Then, the transformed Baldwin-Barth turbulence model given by Equation 
(21-86) is written as follows. The details of transformation and mathematical ma- 
nipulation is provided in Appendix I. 
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where U and V are the contravarient velocity components 
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Furthermore, the expression for S in the last term of Equation (21-92) is 
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4 2b; 
where the coefficients a;, a2, a3, b1, be, bs, C1, C2, cs and c4 are given by (11-2102), (11- 
210b), (11-210c), (11-2112), (11-211b), (11-2116), (11-2122), (11-212b), (11-212c), 
and (11-212d), respectively. The transformation of Equation (21-87) is carried out 
in a slightly different fashion, the details of which are given in Appendix I. The 
transformed equation is given as 
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In order to simplify the development of finite difference equations, Equation 


(21-99) is expressed as 
OR 
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The term M, which represents the convection and diffusion of turbulence can be 
written as 
= 3 
M = M’ + MẸ + M; + MẸ + M} (21-105) 
where 
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The initial and boundary conditions are set according to: 


(1) Initial condition: A value of Rer = (Rer)o < 1 is specified within the domain 
to initiate the solution. 


(2) Boundary conditions: The value of Rer is set equal to zero at the solid surface, 
and it is set to its initial value of (Rer)o at the inflow. At the outflow, 
extrapolation from the interior domain is used to specify Rer at the boundary. 


21.4.2.2 Spalart-Allmaras One-Equation Turbulence Model: The 
Spalart-Allmaras model solves a transport equation for a working variable v which is 
related to the eddy viscosity. The governing equation is derived by using empiricism, 
dimensional analysis, Galilean invariance, and selected dependence on the molecular 
viscosity [21-12]. The transport equation is expressed as 


T LIV (war) + cule)" + ansa - fa) 





Chi v? 
E [efe - 2 Iz + fa(Aq)? (21-111) 
where the eddy viscosity is given by 
Wo Vfa (21-112) 
and 
fu x3 x z (21-113) 


(21-114) 


x 
I 
Tis 
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Various functions and constants appearing in Equation (21-111) are defined as 


v 
g^ (21-115) 


where d is the distance to the wall, « is the von Karman constant, S is the magnitude 
of the vorticity, 








ôv Au 
5 =| - 2 (21-116) 
and " 
v= 1 —- .—c—- 21-117 
f 3 1 + XÍa ( ) 


It should be noted that expression (21-66) has also been used for S. The function 
fu is 





l4 cfs i 
(r) = 21-118 
40 = c (555 (21-118) 
. where 
g =r + Clr? — r) (21-119) 
and i 
r= Ea (21-120) 


Large values of r should be truncated to a value of about 10. The function fi; is 
given by 

fa = es exp(-cux?) (21-121) 
and the trip function fa is 





2 
fu = €ugir exp E (&) (d + | (21-122) 


The following are used in Equation (21-122): 
d;: The distance from the field point to the trip which is located on the surface. 
wt: The wall vorticity at the trip. 
Aq: The difference between the velocities at the field point and trip. 
9: 9 = min[L0, Ag/u, Az], where Az is the grid spacing along the wall at the 
trip. 

The constants used in the equations above are 


o= Cy = 0.1355 Cio = 0.622 


C3] to 


enn = + (1+ c)/o Cy2 = 0.3 €u3 = 2 k = 0.41 


Cy = 7.1 Ci = 1.0 C2 = 2.0 Ci3 = 1.1 Ci = 2.0 
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The Spalart-Allmaras turbulence model given by (21-111) can also be written 
as 





gv (- 1%) Y. vt 007] - (rez)? 


y2 
+ ca (l — fi)vS — [cus fe = 2 fal H + fa(Aq) (21-123) 


The first two terms on the right-hand side of Equation (21-123) can be expanded 
and recombined to provide the following ines 
dv = (= + Coz 
dt — 





) [Viv +) - V9] + Lv n 


+ on(t = fS - [ese - fa] [ET + sat? (21-124) 


Either one of Equations (21-123) or (21-124) can be used. 
The initial condition for 7 is specified to be zero up to a value of 1/10, that is, 
v = 0 — 1/10. The boundary conditions are set according to: 


(1) At the inflow, v = 7, 
(2) At the solid surface 7 = 0, and 


(3) At the outflow, extrapolation is used. 


21.4.2.2.1 Nondimensional form: Using the same nondimensional variables de- 
fined previously, the nondimensionalized Spalart-Allmaras turbulence model is given 
by the following equation. Again, as in the previous model, the notation of asterisk 
has been dropped. 


£ - (EL) C28) vean- (gh) (S) eee 
+ eun(1— fa)Sv — (a) (Cui fu) (5) 


+ (G) (3) (1 — fi) foe + fa] (5 ) + Res fa(Aq)! (21-125) 


which corresponds to Equation (21-123) and 


T > (Rez) CR) wem vm (iL) ev 





katei G) (cur fu) EY 
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KA 


+ (z-) (3) [(1 — fre) fva + fia] EY + Resfu(Aq? (21-126) 


corresponding to Equation (21-124). 


21.4.2.2.2 Spalart-Allmaras turbulence model in computational space: Using 
the expressions for Cartesian derivatives given by (21-80) through (21-91), the 
Spalart-Allmaras turbulence model given by (21-126) in generalized coordinates 
is written as 











Ov, y (1 x. 
al Sort She w 
ar LETS yc mt (1 — fa)Sv - (x) 2 y 
+ 265 =a dean 4 9 5z 8€ 92 an Chi 2)OV Rea CuiJw F 
1 p? 
+ Gu) (B) [a rae] (F) + ne thor ^ Qum) 
where 
ML ta Rud (21-128) 
"BE tho "5E $n 
and, as defined previously, 
a = E+E (21-129) 
bo = utu (21-130) 
cs = Nz t+ Ey (21-131) 
_ E 06. 0€ Oy 
n = 3 +15 "Ju eet" von (21-132) 
One. 0n ,,9 D 
9 = & ae tn Se TÉ wt oly (21-133) 


The transformed Spalart-Allmaras turbulence model given by (21-125) is 
ov jV , y _ 1 1+ Che da ða 98 B 
m'a By si ) lez ^m +age NE) 


ODTÜ KÜTÜPHANESİ 
M, E, T. U. LIBRARY 
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l Y Che OA OA OB 0B 
-( ) e» fe a ua Hage | 





Rew 
1 Col U 2 
* (x) (=) (1 — fre) fo2 + fal (2) = (a) (Cur fw) (5) 
+ en(1 — fa) SP + Reco fa(A^a) (21-134) 
where 
àv Ov ov 
v ov ov 
B = ay E +n 27 (21-136) 
8v àv 
a = (v-V) (ez: t ng] (21-137) 
ov ov 
Equation (21-134) is now written as 
A =M4+P4+D4+T (21-139) 


where 


(a) The term M represents the convection and diffusion of turbulence and, due 
to numerical consideration, is decomposed as follows: 


M = M'+ M? +M’ (21-140) 


Each term in expression (21-140) is defined as 





ov ov 
Es - (vs XU 5.) (21-141) 
— (1) (litte ða, da Ə , OB 
M = (=) ( 5 ) lex 8€ TUR ðn tt 3E + Tv On } (21-142) 
= 1 Cp2 0B 
M? = See —(v van [e thon Mage tn "on j (21-143) 
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The terms M? and M? can be further regrouped as terms including £ deriva- 
tives and terms involving 7 derivatives as follow 


Mi = (a) (=~) («32 a) (21-144) 
«(unm (en) ans 
Mi = - (a) w+) (eoe RES a) (21-146) 
M3 = - (az) 2 S8 (y, 4. y) ("5 = int ee (21-147) 


(b) The second term in Equation (21-139) represents the production of turbulence 
and is 


"U 
| 


= cn(1— fo)Sv + (Cra ijs e — — fiz) foe ET 


Co (1 - fo)v [s T map. = Cn (1 — f2)0S (21-148) 


1 v 
Note that S = S + Ra (=z) fw 


(c) The third term D represents the destruction of turbulence and is given by 


p=-(qis) (we me) mum 


(d) And, finally, the last term in Equation (21-139) represents the trip term which 
causes/ triggers transition from laminar flow to turbulent flow. 


= (Rec) fa(Aq)? (21-150) 


21.4.3 Two-Equation Turbulence Models 


It was pointed out previously that the convection of turbulence is not modeled 
in zero-equation models. Therefore, the physical effect of past history of the flow is 
not included in simple algebraic models. In order to account for this physical effect, 
a transport equation based on the Navier-Stokes equation may be derived. When 
one such equation is employed, it is referred to as a one-equation model which was 
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discussed in the previous section. When two transport equations are used, it is 
known as a two-equation model. 

Complex flowfields which include massively separated flows, unsteadiness, and 
flows involving multiple-length scales occur frequently in fluid mechanics applica- 
tions. In these types of complex flows, the lower order turbulence models, that is, 
zero-, half-, or one-equation models, become very complicated and often ambiguous. 
Two-equation models are developed to better represent the physics of turbulence 
in these types of complex flowfields. In this section, k-e and k-w two-equation 
turbulence models are reviewed. 


21.4.3.1 k-e Two-Equation Turbulence Model: A commonly used two- 
equation turbulence model is the k-e model. The partial differential equations are 
derived for kinetic energy of turbulence (k), and the dissipation of turbulence (e), 
where 

= 3 [u? + v? + w?] (21-151) 


Ou V / Ou; 
«= (24) (2) (21-152) 


Since it is important to identify the terms involved in these equations, and, for 
the benefit of those who are not familiar with the origin and derivation of these 
equations, the details of the kinetic energy of turbulence equation and interpretation 
of terms are provided in Appendix K. Now, the standard k-e two-equation model is 
expressed by the turbulent kinetic equation 


and 


dk ð Ok 
Pu = da; (o A) Dy. d + Pk — pe (21-153) 
and the dissipation rate equation 
de ð mN ôe € e 
— = -— L—-j— Pe — Cop -154 
Pat 02; (+ 3 =| ior d" Veto) 


where P, is the production of turbulence defined as P, = ri; (Qu;/Oz;). 

The turbulent kinetic equation and the dissipation rate equation given by (21- 
153) and (21-154) can be written in an expanded form in Cartesian coordinates for 
two-dimensional problems as 


ok Ok ok\ 8 IAN Ok ð ðk n" 
(aea az" 25) = 2 5) E] s e) a]*^ d 


(21-155) 
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and 
de ud 44%) e 9 (my oe], 2 [fpa my e 
(Frage test) Oe (n+) FE) + (ut) 


€ 2 


€ 
+ caP, kp E (21-156) 
where typical constants are 
Sk — 10 o,— 13 co =144 ca= 1.92 


and the turbulent viscosity is related to e by 


2 
He = peu — (21-157) 


where c, = 0.09. 
Before proceeding further, it is important to note that the left-hand side of the 
k and e equations, given by (21-155) and (21-156), can also be written as 


ô 0, LR 

at (pk)+ On (puk) 4 by (57k) = RHS; (21-158) 
and 

ô Orpa nee 

$i (pe) + an (pue) + 5 (pvc) = RHS: (21-159) 


which is simply obtained by multiplying the continuity equation by k, adding it to 
the left-hand side of (21-155), and subsequently combining the terms. A similar 
approach is used for the e-equation. 

As for the zero-equation models, a two-layer approach is incorporated. That is, 
Equations (21-155) and (21-156) are used in the outer region down to the vicinity 
of the viscous sublayer. For a typical wall flow, that corresponds to a y* range 
of 30-50. Adjacent to the surface, wall functions such as (21-56) are employed. 
The values of k and e must be provided at the interface as the inner boundary 
conditions for Equations (21-155) and (21-156). These values may be obtained by 
the following: if y; denotes the interface between the outer region and the viscous 
sublayer, then 

ke = $ (21-160) 
pch 
where cp is typically 0.164, and 7; is the shear stress at y; location. Similarly, 





i 
Ey-y = cob (21-161) 
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where the length scale may be evaluated by the simple linear relation for the viscous 
sublayer given by 
l2 ky (21-162) 
Once the local values of the kinetic energy of turbulence and dissipation of 
turbulence have been computed, turbulent viscosity given by (21-157) may be de- 
termined. It is noted that correlations such as u?, v?, and w? can be determined 
by semi-empirical relations. For example, 


ul? = 2a;k (21-163) 
v? = 2ask (21-164) 
w = 2(1 — a; — a3)k . (21-165) 


where az and az are structural scales usually assigned values of 0.556 and 0.15, 
respectively. 


21.4.3.1.1 Low Reynolds number k-e model: The difficulty with the standard 
k-e model introduced in the previous section is that the equations become numeri- 
cally unstable when integrated to the wall. To overcome this problem, a two-layer 
approach is implemented, as discussed previously. However, a better approach is 
perhaps to directly integrate the k-e equations through the viscous sublayer all the 
way to the wall. In order to enable the integration to the wall and to improve the 
capability of the standard k-e model, several modifications are introduced. The 
resulting formulation is known as the low-Reynolds number k-e equations. The 
first low-Reynolds number k-e model was developed by Jones and Launder [21-13, 
21-14], and subsequently it has been modified by several investigators. The pri- 
mary modifications introduced by Jones and Launder were to include turbulence 
Reynolds number dependency functions f; and fz in Equation (21-154) and f, in 
relation (21-157). Furthermore, additional terms L, and L, were added to the equa- 
tions to account for the dissipation processes which may not be isotropic. Thus the 
low-Reynolds number k-e equation is written as 


dk ð pu Ok 
d ox: ud Da z 21- 
Pd Oz; | (u p A) z] + Pk — pet Le (21-166) 
and " j " : 
ie, e ut € €. - € : 
Pat x; Ic a) =| +a fi Fey cafa tLe (21-167) 


where the turbulent viscosity is now computed according to 


k? 
He = Phuc (21-168) 
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The following relations are proposed by Jones and Launder. 





A = 10 

fo = 1-0.3exp (—Re}) 
zn | -2.5 
e = 9PHIQORe, 

2 
ovk 
Lr ~ -z (52) 
and 


L 





I 
pas 
D = 
Q 
& tw 
B, g 
— 


Since the early 1970’s several modifications to the low Reynolds number terms 
and the model constants have been introduced. A summary of selected k-e models 
are provided in Table 21.1. 









Model [Ree | A | dz fe 
eae 


225 


Nagano-Hishida 1-0. ame 





[I — exp(—Rer/26.5)] 





2 
Nagano-Hishida vw(1- fa) mm 


Table 21.1 Functions and constants for the k-e turbulence models. 
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21.4.3.1.2 Compressibility correction and axisymmetric consideration: In addi- 
tion to the modifications of the low Reynolds number terms and model constants 
shown in Table 21.1, other modifications have also been introduced into the k-e 
equations, most notably, the compressibiility correction terms and axisymmetric 
considerations. 

The compressibility correction term is typically included by introduction of the 
turbulent Mach number M, defined by 


2k 


a? 


M = 


Now, the k-e equations including the compressibility correction terms are written 


as 

Dama L Gak) = ue) Ok -5 d x 

5i PR) + 5 FUA = 5 IG E) z] + Py B(c + €a) + P/d? + Ly (21-169) 
and 

ô, Q E I4 \ ôe € _é 

à 09 Gy, Pts) = EA I 3 dz; Ca fi hg ca fap tL (21-170) 


The term e, represents the contribution due to compressible dissipation, and the 
term p'd" represents the pressure dilatation term. They are given respectively by 


e = y Me (21-171) 


and 
pd! = —y P, M? +- Y3 pe M? (21-172) 


The following values based on DNS have been suggested. 
10, ^9. — 0.4, and 44-02 


The simplest approach to include axisymmetric correction is to modify the co- 
efficient c, as 
{ 1.44 for planar 2-D flows 
Ca = 


1.60 for axisymmetric 2-D flows 


21.4.3.1.3 Initial and boundary conditions: The initial and boundary conditions 
must be specified for both k and e. The initial condition for the turbulent kinetic 
energy k is specified in terms of the freestream turbulence intensity T;, as follows 


koo = 1.5(T; Uo)? 
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where T; = 0.0001 to 10 percent of mean velocity. The turbulent viscosity is specified 
as 


Hteo = (0.1 — 100) us 


and, therefore, 
k2 


oo 


€oo = PCy—— 
Pts 


The boundary conditions are 
(1) At the inflow, the freestream values are specified, that is, € = €, and k = ky. 


(2) At the solid surface, k = 0, jj, = 0, and € = e, = 0. For some k-e models, a 
value for €y is prescribed. 


(3) At the outflow, extrapolation is used. 


21.4.3.2 k-o Two-Equation Turbulence Model: This two-equation model 
includes one equation for the turbulent kinetic energy k, as developed previously, 
and a second equation for the specific turbulent dissipation rate (or turbulent fre- 
quency) w. 

Recall that the development of the k-equation was based on the governing equa- 
tions of fluid motion. A second approach for the development of a transport equa- 
tion may be considered whereby the equation is constructed based on the known 
physical processes along with dimensional analysis. This is the approach taken to 
establish a transport equation for w. 

The concept of parameter w was introduced by Kolmogorov, which he called 
dissipation per unit turbulence kinetic energy. Recall that the process of dissipation 
takes place at the level of smallest eddies, while the rate of dissipation is the rate of 
transfer of turbulence energy to the smallest eddies. Therefore, this rate is set by 
the properties of the large eddies. Perhaps the simplest physical interpretation of w 
is that it represents the ratio of turbulent dissipation rate to the turbulence mixing 
energy or, alternatively, the rate of dissipation of turbulence per unit energy which 
in essence is an inverse time scale of the large eddies. 

Now, to develop a transport equation for w, several physical processes which are 
observed in fluids are utilized. These processes include unsteadiness, convection, 
diffusion, dissipation, and production.. A combination of these physical processes, 
along with dimensional arguments, yield 


Ow Ow _ 2, 9 Ow 
Pay tp n ppu? + Ba; len z] (21-173) 


where f and c are coefficients to be specified. 
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In comparison to the k-equation, several observations can be made. First, note 
that a production term does not appear in Equation (21-173). Second, a molecular 
diffusion term is also absent from this equation. Inclusion of a molecular diffusion 
term is essential if one intends to integrate the equation through the viscous sublayer 
to the wall. Equation (21-173) is considered as the basic w-equation upon which 
modifications are introduced. 

As in the case of the k-e model, there are several versions of the k-w model. 
Perhaps one of the best known is the k-w model of Wilcox [21-17], which is provided 
below. 

The turbulent kinetic energy equation is given by 


dk 0 | Ok 
ZEE (Houge + Pa~ Bink (21-174) 
dt = 0x; "Oz; 
and the specific dissipation rate equation is 
dw 0 Ow w 
jn Bathe = Ds 2 21. 
P aA EE + oy Pr - Bow (21-175) 
The eddy viscosity is determined from 
Et (21-176) 
=p m 
and the auxiliary relations are 
e = f'uk (21-177) 
ki 
= PA (21-178) 
The constants used in Equations (21-174) and (21-175) are 
5 3 . 9 1 AE 
Weg deu: Pg Emus gnum 


21.4.8.8 Combined k-e/k-w Two-equation Turbulence Model: The 
k-w turbulence model performs well and, in fact, is superior to the k-e model within 
the laminar sublayer and the logarithmic region of the boundary layer. However, 
the k-w model has been shown to be influenced strongly by the specification of 
freestream value of w outside the boundary layer. Therefore, the k-w model does 
not appear to be an ideal model for applications in the wake region of the boundary 
layer. On the other hand, the k-e model behaves superior to that of the k-w model 
in the outer portion and wake regions of the boundary layer, but inferior in the 
inner region of the boundary layer. 
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To include the best features of each model, Menter [21-18] has combined different 
elements of the k-e and k-w models to form a new two-equation turbulence model. 
This model incorporates the k-w model for the inner region of the boundary layer, 
and it switches to the k-e model for the outer and wake regions of the boundary 
layer. Two versions of the model introduced by Menter are referred to as the baseline 
(BSL) model and, a modified version of the BSL model, the shear-stress transport 
(SST) model. It has been shown that the SST model performs well in the prediction 
of flows with adverse pressure gradient. 

The baseline model is identical to the k-w model given by Equations (21-174) 
and (21-175) in the sublayer and the log layer of the boundary layer, and it gradually 
switches to the k-e model given by Equation (21-153) and (21-154) in the outer wake 
region. These equations are repeated here for convenience. The k-w equations are 


8, ð RUNE: Ok) .- S 
5, 09 + Jz P% k) = öz, [o + Oki Ht) Ed + Py — B'pko (21-179) 
and 
9, Ou s ox ð Ow w 23 
a; 0 + Ba, P w) = dn, lo + 24114) =] + ory P iw (21-180) 
where k 
Me = PT (21-181) 
and 
e— B'uk (21-182) 
Furthermore, the following are used 
0k —0', 01-70, =b, and œssa 


The k-e equations are 


DTE IN Ok - 
a, 9 WE 8z; PF) = Bz; | (u+ A) =l + Py pe (21-183) 
and 
0, ne _ ð &) dc € .e 
5, 79 + J Pm e) = Bz; I = =| + CaF P; — cap k (21-184) 
where 
k? 


He = peu (21-185) 


o ARM —X—À B—— E Chapter 21 


and : 

pr: (21-186) 
Now, in order to combine the two sets of equations, the k-e set is transformed into 
a k-w formulation using relation (21-182). In the process of this transformation, 
two additional terms appear in the new w-equation. One of these terms is a cross- 
diffusion term and the other occurs if o, and c, are not equal. It has been shown 
that the second term is relatively small and does not affect significantly the solution. 
Therefore it is neglected. Subsequently, the transformed k-e equations written in a 
k-w formulation are eas as 


à x. 
2 (gk) + 2 2 (pa; ve |o toan) 2. =] + Pe- p'pku (21-187) 
and 


Fe) + 5 Pyu) = x [u+ aama] 


1 ðk ðw 


Yp, — 2p ow 
+ a2 — Py — Bo po? + pa, DIOS 


k 


The relations between the coefficients in the original k-e equations and the trans- 
formed set are established as follows 


(21-188) 


1 
On = n gua = E f = B'(ca — 1) 
BY =c, Q2 = (ca — 1) 


Now the two sets of equations are combined by the introduction of a blending 
function F. The k-w set given by Equations (21-179) and (21-180) is multiplied 
by F, and the transofmred k-e set given by Equations (21-187) and (21-188) is 
multiplied by (1— F). They are subsequently added together. The blending function 
F is designed such that it will be equal to 1 in the vicinity of the wall region, thereby 
activating the k-w model, and it will be zero away from the wall, thus activating 
the transformed k-e model. 

The resulting combined k-w/k-€ EE model is given by 


8 
bont ZOU =i [ut omige|+ P-e nas) 
and 
Boss ce PRU. à 
ài (pw) + x; (Piw) = Ox; le 2-4 + 


1 Ok dw Du. aU ee 


HER 2 j0z; k 


(21-190) 
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where the production of turbulence as defined previously is given by 


Ou; 


P; = T; — 
"On; 


(21-191) 





Ox; 9 


The constants appearing in Equations (21-189) and (21-290) are expressed in a 
general compact form by 


Ou, ô 2.0 2 
ny = me (5 is + A Sis 54) - gos 


$ — Fiji + (1- F1) (21-192) 


where $; represents the constants associated with the k-w model (when F; = 1), 
and $s represents the constants associated with the k-e model (when F; = 0). 

The difference between the baseline model and the shear stress transport model 
is in the definition of turbulent viscosity and the specification of constants. 


21.4.3.3.1 Baseline model: The turbulent viscosity for the baseline model is de- 
fined by 


i p. (21-193) 


and the constants for $ are set as follows. 
The values of the constants for set ¢, are specified according to 


Oki = 0.5 , Owl = 0.5 , i= 0.075 
B' =0.09 , k = 0.41 


a= a TETTE 

M the constants for $ are 
1 

On = — — 100 , Ou. = z = 0.856 , 8, = 0.0828 
€ 


B' —0.09 , k= 0.41 


In addition, the following definitions are used. 


F, = Tanh (arg) (21-194) 











Vk e) d (21-195) 


arg, — min [ax ( 0099y ' wy)? CD? 
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where y is the distance to the nearest surface, and C D,, is the positive portion of 
the cross-diffusion term 


i . 
C Dy, = max |2p0,5— Oke ge 1079 


Note that far away from solid surface as y — large, term arg; approaches zero 
due to 1/y and 1/y?. Furthermore, the three arguments within arg, represent the 
following: 


1) The term (Vk/0.09wy) is the turbulence length scale divided by y. It is 
equal to 2.5 in the log layer of the boundary layer and approaches zero at the 
boundary-layer edge. 


2 


— 


The term (500v/wy?) enforces F; to be one in the sublayer region and 1/y in 
the log layer. Recall that w is proportional to 1/3? near the surface, and it is 
proportional to 1/y in the log region. Therefore, 1/wy* will be constant near 
the surface and will approach zero in the log layer. 


3 


— 


The third term (4po,sk/C Dy?) is included to prevent solution dependency 
on the freestream. 1t can be shown that, as the boundary-layer edge is ap- 
proached, arg; as well as F become zero and, therefore, the k-e model is 
utilized in this region. 


Transition to turbulence is activated by phasing in the production term, typically 
over a few grid points at a specified transition point. Note that the production term 
is turned off in the laminar portion of the flow. 


21.4.3.3.2 Initial and boundary conditions: The following freestream conditions 
are recommended. 


u 

Woo = mọ ; Vin = 10" , koo = Vt Uoc 
where 

l<m<il0, 2<on<5 , 


and L is a characteristic length of the problem, for example, the length of compu- 
tational domain. At the solid surface, the boundary condition for w is set according 


i 6v 


o= Ai? 
where Ay; is the distance of the first point away from the wall. This boundary 
condition is specified for a smooth wall, and Ayt must be less than about 3. At 
the outflow boundary, extrapolation is used. 


(21-197) 
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21.4.3.8.8 Shearstress transport model: The turbulent viscosity of the BSL 
model is modified based on the assumption that the turbulent shear stress is pro- 
portional to the turbulent kinetic energy in the logarithmic and wake regions of the 
turbulent boundary layer. The amount of computed turbulent viscosity is limited 
in order to satisfy the proportionality requirement. This limitation on the turbu- 
lent viscosity is imposed on by the turbulent kinetic energy, and thus the turbulent 


viscosity is defined as 
ak 


= max(ayw , QF) 
where a, = 0.31, Q is the absolute value of vorticity Q = |8v/0z — 0u/0y|, and 
F, = Tanh (arg?) with 


7 (21-198) 


Vk d (21-199) 


0.09wy’? wy? 
The constants for ¢; and @¢2 are identical to the baseline model, except cj, 
where 


arg, = max 2 


Or = 0.85 


21.4.3.3.4 Compressibility correction: Similar to that of the k-e model, the com- 
pressibility correction term is included by introduction of turbulent Mach number. 
Thus, Equations (17-189) and (17-190) are formulated as 


POTRETE ok 
a 09 + au PYA = $z [utom gE] +r 


- B*pwk [1+ aM?(1- F))] + (1 — Fi)" (21-200) 
and 


TENERE: ay 
HO) x09) = g [urou 2] 


"Bm Ra peg oe” 
ter Py — Bpw* + 2(1 Fi)P ou ~ dz, Oz, 


+(1-F)6'a, M}pu?~—(1— Fp" (21-201) 
t 
where the blending function F; is defined as in the original model. 
21.4.8.8.5 Nondimensional form: The nondimensional variables defined previ- 


ously are used to express Equations (21-189) and (21-190) in a nondimensional form. 
Again, the notation of asterisk has been dropped. However, the appearance of the 


—————————— SAI 222 8633 


Reynolds number in the equation is a good indication of nondimensional equation. 
The nondimensional k-c/k«» two-equation model is written as 


dk (5 Ok S) = 1 8 


^p utu + (u+ ont) E 
P Oe ar By) Rea dn |V T I Bg 
1 2 ðk E 
+ Ren ay |o d -B- pwk (21-202) 


dy C E A 
PH =? ôt “On "Oy ^ Res OX Bord die Ox 





1 


à Ow 
+ Re ay lu + Ou wd | + 2p(1— Fijos — 3i 


Ok dw " Ok Ow 
Ox Oz | Oy Oy 


4 or P. - Bp? (21-203) 


where the production term P, is determined from 


| u Ou ðu; 2 7 Our Ou; 
P; = Se, = p (Fe + an 3 ou mt) - 3 jou t On; (21-204) 


21.4.3.3.6 k-e/k-w Turbulence model in computational space: Using the expres- 
sions for the Cartesian derivatives given by (21-88) and (21-89), the k-equation in 
the computational space becomes 


Ok Ok yok 1 ð ð 

Pay = ~PU Gg ~ ° 24 + Rents 5g (MUS) (KX) + ne (MUS) (XO) + 

be zg (MUS) (KY)] +n 2 a; (MUS) (KY) je Pr- Souk (21-205) 
where 

MUS = uou (21-206) 

Ok Ok 
KX = & ge +5, (21-207) 
KY se ee (21-208) 
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Similarly, the w-equation is given by 


= —pU—-—pV—+ 


Ow dw Ow 1 { 
Pi 8€ an Re, 


à 
& gc (MUR) (020] 

2 MUR)(OX d MUR)(OY g MUR)(OY 
+e gy (MUR) (OX)] + & 5c (MUR) (OY)] + ny 5 (MUR) (OY) ) 


+ [20 — Fes =] (KX) (OX) + (KY) (OY)] — Bow 


(21-209) 
where 
MUR = uou (21-210) 
ô 
OX =& o +1: » (21-211) 
w Ow 
OY = & 8€ + Ty on (21-212) 


The production term Pj given by (21-204) is first expanded and subsequently trans- 
formed to the computational space as follows. 


i du ô 0 à 
SEE = ~ Teg, + Tey (S us 5) +w ay oa 
where 
1 4Q0u 2 ov _2 
1 Qu Ov 
ey oa ee 21-215 
Try (a) di (5 m z) ( 
1 40v 20v 2 
m dh. SOUNA 2 1-21 
Tyy (az) p (5 Oy 3 z) - 3^* tale) 


Utilizing the expressions (21-88) and (21-89), the shear stresses defined by (21-214) 
through (21-216) become 


1 Ou Ov 2 
Tr, = (zz) ls (e: 8€ +1652) - 3 (Gee + Ty n =) ee 3^* = Teg (21-217) 
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i Ou Ov 
_ - 21-21 
T (zz ) me lege Tou tex i 5. Tin ae 


1 4 dv Ov Qu 2 
tw = (Fez) $ (ez tna) - (eat 5) o d 


Now, the production term becomes 


Py = Taza 


Ov ) 
+ beget ge JE (ase ear) 


ðu 
(Ez Tee + & ren) oe p (nz Tee + Ny Ten) n 


ô Ov 
+ (& Ten + £y Tn) o + (nz Ten + Ny Tm) ôn 


E (P1)? = EF (Pa) Ê z x (P3) 2 ps (Ps (21-220) 
where 
Pl = £ Tee + Tm (21-221) 
P2 = Na Teg + Ny Ten (21-222) 
P3 = & Ten + &y Ton (21-223) 
P4 = moe + Ny Tm (21-224) 


21.5 Numerical Considerations 


Numerical simulation of turbulent flow may be divided into three categories. 
The simplest and most practical approach in use at present is the solution of the 
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time-averaged Navier-Stokes equations along with a turbulence model. The second 
level of sophistication includes the simulation of time-dependent, large-scale eddy 
motion. In this approach, the effects of smaller scales are included by turbulence 
models. The third category includes the direct numerical simulation of all important 
scales of turbulence. 

The second and third categories are presently in the infancy stage. One of the 
difficulties associated with these methods is the limitation of present day comput- 
ers. However, with the advancement in computer technology, large eddy and full 
turbulent simulations seem possible for practical problems within the next decade 
Or SO. 

Currently, the system of equations such as (21-27) along with some turbulence 
models are solved for engineering applications. The solution schemes for the system 
of equations given by (21-27) have been discussed in the previous chapters. One 
` important point to emphasize here is grid resolution. In order to accurately resolve 
the feature of turbulent flow near the surface, a grid point within the laminar 
sublayer must be included. Generally, this corresponds to a grid point within y* of 
2.0. 

The inclusion of a turbulence model into the system of governing equations 
may be accomplished either explicitly or implicitly. One may solve the governing 
equations [such as (21-27)] implicitly and use the turbulence quantities from the 
previous time-level, i.e., explicitly. On the other hand, the governing equations, 
along with turbulence equations such as (21-155) and (21-156), may be solved si- 
multaneously, i.e., implicitly. Obviously, a fully implicit formulation will increase 
the size of the Jacobian matrices and will require major modification of the (lam- 
inar) Navier-Stokes code. On the other hand, explicit treatment of turbulence is 
simple and straightforward, and modification of an existing Navier-Stokes code to 
include turbulence can be accomplished with ease. 


21.6 Finite Difference Formulations 


The one-equation and two-equation turbulence models in different forms given 
in Sections 21.4.2 and 21.4.3 are parabolic equations. Therefore, various numerical 
methods discussed previously for the solution of parabolic equations can be used to 
develop solution algorithms. 

To illustrate the development of the numerical schemes for the turbulence mod- 
els, both explicit and implicit formulations are considered. However, the implicit 
formulations will be limited to the one-equation models. The turbulence models 
in nondimensional form and in the generalized curvilinear coordinate system are 
considered to provide the FDEs for general applications. 


70 Chapter 21 


21.6.1 Baldwin-Barth One-Equation Turbulence Model 


A first-order upwind approximation is used for the convective terms, and a 
second-order central difference approximation is used for the diffusion terms. The 
finite difference expression for each term is developed separately due to the com- 
plexity of the equation. The first-order upwind scheme can be easily extended to a 
second-order upwind scheme. 


21.6.1.1 Implicit Formulation: The Baldwin-Barth one-equation turbu- 
lence model written as aR 
a MP (21-225) 
is used to illustrate the development of an implicit scheme. Equation (21-225) is 
written as agn 
(F) = M"! + prt} (21-226) 
ôt 
Since the governing equation is nonlinear, a linearization procedure similar to 
that introduced previously, for example, by relation (11-81), is used to provide the 


following. "T 


Mr =M" + zp AR=M"+M"AR (21-227) 
and ap 
P! = pra JK AR = P” + PAR (21-228) 
Now Equation (21-226) is written as 
SE — (M + P)AR = (M? + P") (21-229) 


Recall that the term M can be decomposed into terms involving € or ņ derivatives 
as defined by (21-107) through (21-110). Thus, M is expressed as M = M; + Ms. 
Now, Equation (21-229) is written as 
AR VA Eva D n n n 
sy (M; + M, + P) AR = (Mẹ + M} + P”) (21-230) 
or 
[1 - (Me + M, + P)At] AR = (MẸ + Mọ + P")At (21-231) 
and decomposed as 
[1 — AtM,| |1- At(M, + P)| AR = (M? + My  PP)AL (21-232) 
Equation (21-232) is solved sequentially be the following equations. 
(1— At M,)AR* = (MẸ + Mj + P")At (21-233) 


and 
[1 - AM, + P)| AR = AR* (21-234) 
The finite difference expressions for each one of the terms appearing in Equations 
(21-233) and (21-234) are developed separately. 


The convective term M! given by (21-106) is LAppredAe by a first-order 
upwind formulation as follows. 


bu upeH QR 
M = Use E 
1 fs Rx. 
= - [5e + Iai) [^ eng) 


1 auus deu 
+ gU — |U,;l) (^u Za) 
tiv qyp (Ric Baa 
+ 3s Ip (Pe Bus) 
1 
+ 5(Vis - IVul) (Bas Su 22] (21-235) 


If a second-order approximation for the convective term is desired, then the formu- 
lation is as follows. 


lu — —9 en — 
MSU V o 
m dd. n (8g — 4Ri-i + Rio 
- imu nup (Mest Res 
lo. apu [ZE + Anas 7 3Rig 
+ gt - sb ( S ) 


ls o (SR - 4Rij-1 + Ris 
ego e up (Wat Ren) 


las waírBuet4ABga-3R : 
+ 3(Vu Vial) [Sev tee fu An )} (21-236) 


Since all the terms within M? given by (21-107) or (21-108) are similar, the de- 
velopment of the finite difference expression will be illustrated in detail for one 
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term only. Subsequently, the conclusion will be extended to the remaining terms. 
Now, consider the term &,(0a/0£) which is approximated by a second-order central 
difference expression as follows. 


eo bes (21-237) 
Furthermore, recall that from (21-84) 
OR OR 
= lr 1- 
a= nh (ez BE +n a) (21-238) 


The finite difference expression for o can be written as 


Rinj — Ri Roa Raa 
Oed = V l&. ( +j — AM Vni c et (21-239) 


AE An 
and 
reda Beini-BRQaisg 
Oi pg = ns le... (Stew) m Mey ebd] (21-240) 


Substitution of (21-239) and (21-240) into (21-237) yields 


da — 1 Rug Ray Bay -Haag-l 
EzE Sa Ezis ac bas le. (ie) + LZ (T 


Rij — Ri-ij Bag 7 Ray 
CA gs leu (Su e) Ys 4 (rii ad j (21-241) 


A similar expression is developed for £,(808/8£) and is combined with &,(80//80£€) to 
provide an approximation for M? as follows. 


gc cde pee [4 E 
M? = i cum E BE + a) (21-242) 
where 
Ga OF n Pig Thy 
t 8€ + yz 8E = Mays (e. ens * Sys Enaga) E 


coe g &, dict Sus Enga) Rcx zal 
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Bp Hu 


Rp eB. 
Tigy (e. Tei us F Sus M45) s d | 


(21-243) 
Similarly for MŽ, 


CA fuic 


~My e y E yh, y) eS t ii 


Bas = Flags) 


TORT Ae we 
"d d TL TW og taa (21-244) 


AgAn 
The terms in M? are also similar to that of M?. Recall that 
2 
3 M 
M = ml’+ 2) («2 +657) (21-245) 
2 v ðA ðB 
ENS A uiis 25 J 
M, = Re. (v+ i) (n. ôn + Ny 5) (21 246) 
and 
ô 
A= e +n A (21-247) 
OR 
B = Guetta (21-248) 


Now the finite difference approximation for M? following (21-243) can be ex- 
pressed as 


2 YE 


Ry - Rag 
3. (Ez fey + bus Enga) eee 


Mo Chapter 21 


Ras — Rigg 
F Ges Megs + Ens Mrs) Esc us j 


Boites 
m (£s, Tis, as t £y Tyga) [F bJ- ii (21-249) 





AgAn 
Similarly, following (21-244) 


3. 2. V Rij+i — Rij 
ms ee pec (ma tt | (An)? 


R; 
(Mez + hysly a) E BAD | 


+ (Massas + Thes.) 


Bags 7 Raus 
AgAn 


pla. — Ri; see 
m (Masz + Tyson, 4) EI (21-250) 

At this point, the finite difference approximation of the terms in Equations (21- 
233) and (21-234) have been identified. However, the bar quantities such as M 
require some deliberation. 

The procedure for obtaining these quantities is illustrated for M, in detail, and 
the remaining terms can be obtained in a similar fashion. 

Recall that M? is approximated by (21-235) when first-order approximation is 
used. Thus, 


OR 
Mj = ae -5 (Us +t) (58 xU) 


1 Risij — Rig \ 
-(U,; — |U,; ioci HAE S A 21-251 
+ gts - Ih (Bey Be) (21-251) 

Note that, as before, when second-order approximation is used, that is, Equation 
(21-236), the terms at grid points i — 2 and i + 2 are taken to the right-hand side 
of the equation and treated explicitly. Now from Equation (21-251) expressions for 
M, are developed as follow. 








=, OM: 1 1 
Ty um [a0 Ml (ae) 
UE OM; 1 1 
M; TREE SR = ess mun (az) 
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a US - Ul) (z) 
= - [ua (a) 


os aM: B 1 1 
M P iin (30s - |U;;l) (z)] 


Elisis OR 
Now, with all the terms in Equations (21-233) and (21-234) identified, the two 
tridiagonal systems of equations are solved sequentially by any tridiagonal solver, 
for example, by the scheme introduced in Appendix B. 











21.6.1.2 Explicit Formulation: The Baldwin-Barth turbulence model given 
by (21-99) is considered to develop an explicit scheme. Again, as in the implicit for- 
mulation, an upwind approximation is used for the convective term which could be 
either first-order or second-order, and a central difference approximation of second- 
order is used for the diffusion terms. In the following formulation, a first-order 
upwind scheme is used. However, as seen in the previous section, extension to 
second-order is simple and straightforward. 


Ag 


1 
+ zU — |U,l) (Bg zu) 


uu 
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where 


2 2 
Uitij Uil Ui j+1 — Uij-1 
Sy = ay (EEN) en Garces 
Uigij — Ui-lj 
vo, (EH) en, (ES 
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2c, | ———— —- 2 —J M 
s, erige) (P asa, (Sn enn) (MH) 


2 
Vij+l — F tet) 


Ui+ij — Wi-j Vit yj — Vi-1g Uijt1 — Wig-1 Uijg1 — Uij-1 
tana (een (unita — ure 
Eur = AT 1 Viel, — Vi-1y 
2AE 


(21-253) 


4 2¢s,, E (eee = Uj-1 


21.6.2 Spalart-Allmaras One-Equation Turbulence Model 


As in the previous section, an upwind approximation is used for the convective 
terms, and central difference approximation is used for the viscous terms. Since 
the approach is similar to that of the previous section, the final form of the finite 
difference expression and/or equations are provided in the following subsections. 


21.6.2.1 Implicit Formulation: The convective term can be approximated 
by either a first-order or a second-order upwind formulation. The first-order ap- 
proximation yields 
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Ov Ov 
-(vi «vi 
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1 Taiga P 
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diy. jy p ( Pi — n 
+ Hs - up (Passt) (21-254) 


The second-order approximation provides 


— AV L; Via 

|! n J j 

Mi E -5 (Us + n (s at Ex) 
= 1 e X: —Vjiajt4Viij — 3Vig 
m+ FU, - Wap (Dent meu Tu 


1 39; Av it Tij- 
+ 104 Mul (Be gee) 


2An 


Ti jt + Wijti — 3Vig 


1 — 
+ 3(u - Vl) EL (21-285) 


The terms M? and M? given by (21-144) through (21-147) are approximated fol- 
lowing the procedure illustrated in the previous section. 


M = (a) (x2) G +6) 
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T (vU TE + Tsw 2 E +h mee (21-259) 
The production term in this model is simulated by 
P = c Sv(1 = fo) 


where cs, = 0.1355 is a constant, and 
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The destruction term is given by 
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The trip term is given by 
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where 
wt 2 2p 
fu = cugeexp COCA ae (a T 9i dj) (21-276) 
and 
g: = min (on, | (21-277) 


Aq: Is the difference between the velocities at the field point and at the trip point 
at the wall. Since the velocity at the wall is zero (no suction), Aq is simply 
the velocity at the field point. 


wi: Is the wall vorticity at the trip. 
Az Is the grid spacing along the wall at the trip. 
di: Is the distance from the field point to the trip point. 
d: Is the nearest distance from the field point to the wall. 


At this point, all the finite difference expressions are in place to develop a nu- 
merical scheme. The ADI formulation, which was introduced in Chapter 3, will be 
used for this purpose. Recall that the governing equation given by (21-139) is 


S S M+P-D+T (21-278) 
Furthermore, Equation (21-278) is nonlinear due to the terms in M, P, and D. 
Thus, before a numerical scheme is considered, Equation (21-278) must be lin- 
earized. The linearization procedure developed in Chapter 6, and similarly devel- 
oped in Chapter 11, and given by Relation (11-81), is used to linearize Equation 
(21-278). The general expression is 


E" = E" + AAu (21-279) 


where 
_ ôE 
^ ðu 


Now, to develop an implicit scheme, Equation (21-278) is applied at time level n+ 1. 
Therefore, 


A 


ðv n+l 
(z) = M™1 4 pP™ pe a qm (21-280) 


Turbulent Flow and Turbulence Models 83 


where the trip term is treated as a source term evaluated at time level n. Now the 
linearization (21-279) provides 


M" = M". OM Ap = M"+MAv (21-281) 

P! = P+ OP Ay = P" + PAv (21-282) 
= 35A" 7 7 

D" = D*+ oO av = D" + DA? (21-283) 


Now Equation (21-280) can be written as 
Av TH 5 Bat = n n n 
Ay ~ (M+ P- D)Av=M +P"-D°4+T (21-284) 


or 
[1 - M + P — D)At| Av = (M? + P^ — D^ + T')At (21-285) 


Recall that, in the development of finite difference expressions for M, it was decom- 
posed as M; and M,. Therefore, the Equation (21-285) can be written as 


[1 - (M; + M, + P — D)At] Av = (M^ + P^ — D" + T')At 
Now this equation is factored as 


[1 - AtM;] [1 - Ac (M, + P - D)] Av = 


(MẸ + Mọ + P^ - D^ +T”) At 
and solved sequentially with the following equations 
(1— At M) Ao* = RHS (21-286) 
[1 - At (M, +P- D)| Av = Av (21-287) 


As seen previously, each one of these equations, once applied to all of the grid 
points, will result in a tridiagonal system. For example, Equation (21-286) will yield 


(1 — At Mg) JAT} 


aj + (1— ACMQUAN 


The finite difference expressions for the terms in RHS;; of Equation (21-288) have 
already been identified. The “bar” quantities are determined in a similar fashion 
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as that of Section 21.6.1.1. For example, the results for M; are as follows 


ea -5U - lUi!) (x) 


d, = - boo (ae) 


1 1 
aU + sl) (3) 








17! 
M; 


iJ 





i-lj 


Similar expressions can be obtained for the remaining terms. 

Once the tridiagonal system of equations given by (21-286) is solved, the right- 
hand side of (21-287) is known, and, subsequently, it can be solved. The system of 
tridiagonal equations given by (21-286) or (21-287) can be solved by any tridiagonal 
solver such as the scheme discussed in Appendix B. 


21.6.2.2 Explicit Formulation: The explicit formulation of Equation (21- 
127) is similar to that of (21-252) and is given by 
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21.6.3 Two-Equation Turbulence Models 


Similar to the one-equation models, the convective terms are approximated by 
either a first-order or a second-order upwind scheme. Second-order central differ- 
ence approximation is used for the remaining terms either at the grid points or at 
midpoints. If central difference approximation is applied at the points next to the 
boundaries, a difficulty will be encountered due to the appearance of points in the 
finite difference equation which are outside the domain. At these points, one-sided 
approximation can be used to eliminate the difficulty. When approximation at the 
midpoints is used, this difficulty does not occur. Both formulations are provided in 
this section. In the following formulations, the first-order upwind scheme is used. 
Furthermore, since the many terms in the k-equation, e-equation, and w-equation 
are similar, only the finite difference equation for the k-equation is provided in 
detail. Subsequently, the extra terms in the e-equation and the w-equation are 
investigated. The k-equation given by (21-205) is divided by p and written as 
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Now an explicit formulation for Equation (21-290) is written as follows 
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where the production term PX, as given by (21-220) is computed as follows 
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To illustrate the midpoint approximation, consider the first four terms on the 
right-hand side of Equation (21-205). Keep in mind that the convective terms 
are approximated by upwind scheme, as before. The four terms of (21-205) under 
consideration are 
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Now consider the finite difference expression for (€,0A/0£) where central differ- 
ence approximation is used. 
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Finally, the finite difference equation for the k-equation becomes, 
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The w-equation has similar terms as the k-equation and, therefore, the FDE just 

investigated can be used for the w-equation by replacing k with w. However, recall 

that the w-equation has an extra term which will be investigated at this point. 
The seventh term of the w-equation given by (21-209) is 
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The finite difference for the expression above is written as 
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21.7 Applications 


Several turbulence models have been introduced in this chapter. The appli- 
cations of these models to three different flowfields are illustrated in this section. 
Computations are performed with several turbulence models for each domain such 
that the solutions can be compared to each other and the experimental data. 

The flowfield is solved either by the first-order and/or second-order flux-vector 
splitting schemes described in Section 14.4.1.2 and given by Equations (14-21) and 
(14-22) or by the modified fourth-order Runge-Kutta scheme described in Section 
14.4.1.3. 


21.7.1 Shock/Boundary Layer Interaction 


The oblique shock generated at a compression corner and its interaction with 
turbulent boundary layer developed near the surface occurs in several engineering 
applications. Thus, the prediction of such a complex flowfield is essential for design 
and analysis purposes. 

Consider now a two-dimensional compression corner in a supersonic flow. Sev- 
eral physical observations have been made with regard to this type of flowfield as 
illustrated in Figure (21-7). (1) At the corner a shock will form causing rise in pres- 
sure. From basic compressible fluid mechanics it is well known that the strength of 
the shock and subsequent rise in pressure will increase as the angle of compression 
corner, and/or the free stream Mach number is increased. (2) A consequence of the 
flow separation is the appearance of a stagnation point at the reattachment location. 
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Figure 21-7. Illustration of a typical supersonic flow at a compression 
corner. 





Note that at this stagnation point higher values of pressure and heat transfer are 
produced. (3) The shock/turbulent boundary layer interaction is highly unsteady. 
The shock typically oscillates back and forth, creating large pressure fluctuations 
at the surface. This physical phenomenon has been shown by several researchers, 
for example, Dolling and his co-workers [21-19 — 21-21]. 

A detailed flowfield survey for several compression corners have been conducted 
by Settles [21-23]. The flowfield corresponding to a compression corner of 24? is 
used in the application presented in this section. The flow conditions are specified 
as follow, Mæ = 2.85(u., = 1875 ft/sec), To = 180 °R, and a constant wall tem- 
perature of Ty = 498 ^R. The Reynolds number base on 6o, where fo = 0.83 inches, 
was 1.33 x 10°. 

The grid system is composed of 101 x 81 grid points as shown in Figure 21- 
8. The grid points clustering near the surface and at the compression corner is 
enforced. The y* for the first grid point off the wall, (i = 2), varied from about 1 
at the inlet to about 4 at the downstream of reattachment point. It is important 
to reiterate that for turbulent flow computations, a y* of less than 5, preferably a 
y* of 1 ~ 2 are necessary. Satisfying this requirement will ensure that at least one 
grid point is located within the viscous sublayer. 
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Figure 21-8. The grid system for the compression corner. 





The initial conditions at the inlet can be obtained by several methods. For ex- 
ample, a computation can be conducted over a flat plate and at a location where 
the computed boundary layer thickness matches the experimental value, the data, 
including the turbulent viscosity distribution, is stored to be used as the inlet con- 
ditions. 

The no slip boundary condition and the constant wall temperature are imposed 
at the surface, that is, the lower surface. At the upper boundary, the free stream 
conditions are specified. Finally, at the downstream boundary, a first-order extrap- 
olation is used. The solution is initiated by imposing the specified flow condition 
at the inlet over the entire domain. 


Comparisons of the computed surface pressure distributions with the experi- 
mental data are shown in Figures (21-9) and (21-10). The solutions shown in these 
figures are obtained using the Spalart-Allmaras one-equation model and the Jones- 
Launder two-equation k-e model, respectively. 


Recall that one of the physical phenomena of the shock/turbulent boundary layer 
interaction identified previously is the unsteadiness of this interaction. However, 
the solutions obtained in this example are developed for steady state. There are a 
couple of reasons for doing so. First, the experimental data used in comparisons 
are provided as their mean values. Second, and perhaps more important, is the 
formulations used in the computations. Since an averaging process is used in the 
development of the Reynolds Average Navier-Stokes equation, the occurrence of 
any high frequency processes within the domain cannot be predicted. In order to 
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Figure 21-9. Surface pressure distri. Figure 21-10. Surface pressure distri- 
butions, Spalart- Allmaras butions, Jones-Launder 
Turbulence model. k-e Turbulence model. 





accurately predict the unsteady nature of this problem, large eddy simulation or 
direct numerical simulation should be used. Nonetheless, the RANS equations com- 
plemented with a turbulence model can predict the overall behavior of the flowfield. 
Comparisons of the velocity profiles at the compression corner and downstream of 
the corner are shown in Figures (21-11) and (21-12). The computed values compare 
reasonably well with the experimental data in light of the complexity of the flow- 
field. Observe that the reverse flow at the corner is predicted fairly well, as shown 
in Figure (21-11). 

It should be noted that, due to the uncertainties in experimentally measured 
turbulence quantities, the approximations used in the RANS equations, and the 
uncertainties in turbulence models, a prediction with a range of +20 percent can 
be considered reasonable for complex flows involving turbulence. 


21.7.2 Two-Dimensional Base Flow 


Near-wake flowfield behind a finite thickness base generated when two streams 
of turbulent flows merge occurs in several applications. The complexity of flowfield 
is manifested if the flow is supersonic. In this type of flowfields, the upper and 
lower streams separate at the corners of the base forming turbulent shear layers. 
These shear layers interact with the expansion waves near the corners and undergo 
recompression and reattachment downstream of the base. Near the base, there is 
a low speed-recirculating region approximately at a constant pressure bounded by 
the shear layers. The base flowfield is illustrated in Figure 21-13. 
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Figure 21-13 Illustration of the two-dimensional supersonic base flow. 
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The objectives of this exercise are to numerically simulate this complex flow- 
field and, in the process, to evaluate the performance of several turbulence models 
discussed previously. 

The flow properties at the upper and lower streams are designated by subscripts 
1 and 2, respectively, and they are specified as 

Mi = 2.56, u, = 584.13 m/sec, pi = 27,566 N/m? 
and 
M = 2.05, uz = 523.83 m/sec, p; = 61,110 N/m2 
The corresponding total pressure and total temperature are 517KP and 299 K, 
respectively. The properties of stream 1 are used to nondimensionalize all the flow 
properties. 

The domain of solution is shown in Figure 21-14, where the origin of the coor- 

dinate system is located at the upper corner, as illustrated. 


Computations are performed on a multi-block grid system where the domain is 
divided into three blocks, as shown in Figure 21-15. The governing equations are 
solved on each block by communicating the solutions at the connection boundaries 
of blocks. Since this approach has not been previously introduced, some deliberation 
is warranted. 

Consider a domain shown in Figure 21-16. For simplicity assume that it is 
required to decompose the domain into two blocks, as shown in Figure 21-17, where 
transfer of data between the two blocks is accomplished through the connection 
boundary. The overlapping lines CD and EF are used as the connection boundary 
between blocks 1 and 2. Note that line EF is identified by i = IM1 for block 1, 
and i = 2 for block 2. Similarly, line CD is identified as i = IM1 — 1 for block 1 
and line i — 1 for block 2. 

The solution procedure consists of two steps at each time level. Solution begins 
in block 1 where, for the first time, the specified initial condition is used to initiate 
the solution. Subsequently, solution proceeds to block 2 where the data along line 
CD (i — IM1 — 1) of the block 1 which has just been completed are used as the 
boundary data along i — 1 of block 2. That is, the boundary values are set according 
to 


()Bleck 2 = (P1-13) Block 1 
(u1)Block 2 = (urmı-15)Block 1 
(My)Block 2 = (vmn-iu)Block 1 


(Cu,)Block 2 = (Ctmmas)Block 1 
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Figure 21-14. Domain of solution for Figure 21-15. Illustration of the multi- 
the base flow. block domain of solu- 
tion. 
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Figure 21-16. Domain of solution. Figure 21-17. Decomposition of the 
domain into two blocks. 





Once the solution for block 2 is completed, the computed values of the variables 
along the line EF of block 2 are used as the boundary values along line EF(i = 
IM1) of block 1, and solution of block 1 is initiated for the next time level. The 
specification of boundary condition from block 2 to block 1 is set according to 
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Figure 21-18. Initial rectangular grid Figure 21-19. Adapted grid system, 
block 1: 44x32, block block 1: 44x70, block 
2: 100x131, block 3: 2: 100x245, block 3: 
44x18. 44x50. 





Computations for the base flow are performed on the rectangular multi-block 
grid system shown in Figure 21-18. This initial computation is performed to obtain 
an estimate of the shear layer in order to generate à more appropriate grid where 
there is a clustering of grid in the shear layer region, and, furthermore, the grid 
lines are more or less aligned along the shear layer. Thus, once a solution on the 
grid system shown in Figure 21-18 is obtained, a second grid based on the flow 
data is generated which is shown in Figure 21-19. Subsequently, the solution from 
the initial rectangular grid is mapped onto the new grid system and used as the 
initial condition. Turbulent flow computations are initiated on this grid. Transition 
to turbulence is specified on the upper and lower surfaces in blocks 1 and 3 at a 
distance of half base thickness downstream of the inflow. Observe that there is 
clustering of grid points along the upper and lower walls and at the inlet location 
in addition to grid clustering along the shear layer. 

As the solution develops, the grid system may be updated in order to better align 
the grid line along the shear layer and to impose more efficiently the grid clustering. 
The inlet flow conditions are generated in a similar approach, as described in the 
previous example. 

Computed streamwise mean velocity profiles obtained by the Baldwin-Barth 
one-equation model, the Spalart-Allmaras one-equation model, and the Baseline k- 
€/ k-» two-equation model with compressibility correction are shown in Figures 21- 
20 through 21-22. The axial locations in Figures 21-20 through 21-22 are at 
xz = 5mm, 25mm, and 55 mm, respectively. Observation of the experimental data 
indicates that the axial locations of z = 5mm and 25mm approximately correspond 
to the separation/recirculation region and the shear layer mixing region, where an 
axial location of z = 55mm corresponds to the recompression/reattachment region. 

The streamwise velocity profiles obtained by the baseline two-equation turbu- 
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lence model with compressibility correction provide the best prediction, as com- 
pared to the experimental data. Comparison of the baseline pressure is shown in 
Figure 21-23. The pressure values obtained from the baseline two-equation model 
with compressibility correction and the Baldwin-Barth model compare well with 
the experimental data, whereas the pressure values predicted by the baseline and 
Spalart-Allmaras models are much lower. 


21.7.3 Axisymmetric, Supersonic, Turbulent Exhaust Flow 


A typical supersonic exhaust flow includes viscous/inviscid interaction, turbulence, 
and possibly it may include thermochemistry. The boundary layers developed in 
the interior and/or exterior walls of the nozzle produce shear layers at the exit 
plane. A complex series of expansion waves and compression waves are generated 
which interact with the shear layer. As a consequence, the so-called barrel shocks 
are formed. Shock reflection from the axis of symmetry may occur in the form of 
either a Mach Reflection or a regular reflection. A typical flowfield is illustrated 
in Figure (21-24). The plume structure is further influenced by the external flow 
conditions, which may be subsonic or supersonic. Other significant physical factors 
which influence the plume flowfield are turbulence and thermochemistry. Typically 
the shear layer is turbulent, and, therefore, appropriate turbulence models must 
be considered in flow simulations. Furthermore, due to large variations in the 
temperature field within the plume, a proper chemistry model should be utilized 
to capture the physics of the problem. In fact, the temperature behind the Mach 
disc could reach as high as the stagnation temperature. Since the primary objective 
in this section is to investigate solutions of several turbulence models without any 
added complexity, the effect of chemistry is not included and the perfect gas model 
is used. Indeed, for the example problem considered, temperatures do not reach 
values for which chemistry consideration would be warranted. 
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Figure 21-24. A typical jet exhaust flowfield. 


Figure 21-25. The grid system for the jet exhaust flow, block 1: 
56x52, block 2: 200x 78 
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Consider now the exhaust flow of an underexpanded nozzle at Mach 2 into a 
quiescent ambient with a pressure ratio of P;/ P, = 1.45, where P; is the jet exit 
pressure and P, is the back pressure. Furthermore, the temperature ratio is given 
as T;/T, = 1.0. The freestream pressure and temperature are 101,325 N/m? and 
293K, respectively. The jet exit radius Rj is 2.54cm. 

A grid system composed of two blocks, as shown in Figure 21-25, is used for 
computations. The grid system includes 56x52 and 200x78 grid points for blocks 
1 and 2, respectively. Due to symmetry, only the upper half of the domain is 
considered in computation. Block 1 extends upstream from the jet exit plane 2.5R;, 
whereas block 2 extends 15R; downstream from the jet exit plane. The outer 
boundary is set at 9R; from the axis of symmetry. The origin of the coordinate 
system is set at the jet exit plane on the axis of symmetry. 

The boundary condition at the jet exit plane is specified, based on the experimen- 
tally measured velocity profile. In the course of investigation, it was concluded that 
the solution is sensitive to specification of boundary conditions set at the freestream. 
In fact, depending on how the boundary conditions were imposed, it was necessary 
in some cases to specify a small external coflow in the order of Mæ = 0.01 or larger 
in order to obtain a stable solution. A detailed investigation of boundary conditions 
for the jet exhaust flow is provided by Hoffmann, et al., [21-23]. 

Turbulent flow computations are initiated by setting onset of transition at three 
grid points downstream of the jet exit plane. The transition is simulated by phasing 
in the production term in the turbulence model over several (e.g., three) grid points. 

The center line pressure distribution obtained with the Baldwin-Barth model 
and a modified Baldwin-Barth model are shown in Figure 21-26, where the numer- 
ical solutions are compared to the experimental data of Seiner and Norum [21-24]. 
The Baldwin-Barth model used in this example corresponds to Equation (21-77), 
with ca = 1.2. The modified version, which includes axisymmetric correction, is 
again Equation (21-77) with simply increasing the value of c4 to 1.47. 

Solutions by Spalart-Allmaras model are shown in Figure 21-27. The Spalart- 
Allmaras one-equation model is given by Equation (21-124) where cp is specified 
as 0.1355. Observe that this model produces an excessive amount of turbulent 
dissipation, and, as a result, the shock cells are dissipated at around z/d of about 
10. The modified Spalart-Allmaras model is the same as given by Equation (21- 
124), except a value of 0.06 is used for cs. Finally, the solutions by the baseline 
k-e/k-w two-equation model are shown in Figure 21-28. Again, as in the case of the 
Spalart-Allmaras model, the original baseline model produces excessive turbulent 
dissipation. Terms are included for compressibility correction (BSL-cc) with a, = 
1.0, ag = 0.4, and a3 = 0.2, and compressibility and axisymmetric corrections (BSL- 
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cca) with o, = 1.0, a2 = 0.4, a3 = 0.2, and o = 0.6. The best solution is obtained 
when both compressibility and axisymmetric correction terms are included. 
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Figure 21-26. Comparison of the pressure distributions along the jet 
centerline, Baldwin-Barth model. 


* Experiment 
--- Original S-A 
| -- Modified S-A 


REUS 
T 
Ee 


EL 


ATI EM. 


= 


rh 


= 
: i ADERAT 
: caren 
wo 


D 





St 
5 


x/d 


Figure 21-27. Comparison of the pressure distributions along the jet 
centerline, Spalart-Allmaras model. 
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Figure 21-28. Comparison of pressure distributions along the jet cen- 
terline, «-e baseline model. 





21.8 Favre-Averaged Navier-Stokes Equations 


Previously in Section 21.3.1, the Reynolds-averaged Navier-Stokes equations, 
based on time averaging, were developed. One of the main assumptions imple- 
mented in the development of RANS is the Morkovin's hypothesis which allows the 
elimination of terms involving density fluctuations. Recall that Morkovin's hypothe- 
sis states that the structure of turbulence is negligibly affected by compressibility up 
to Mach number of about 5 and remains similar to that of incompressible flow. To 
remove this assumption and to increase the range of applications to high speed com- 
pressible flows, the density fluctuating correlations such as p'u’, pv’ are retained. 
The resulting equations are known as the Favre-averaged Navier-Stokes (FANS) 
equations which are also referred to as the mass-averaged Navier-Stokes equations. 
Thus, the terms Favre-averaged and mass-averaged are used interchangeably. The 
objective of this section is to review the development of FANS equations. 

A mass-(or Favre-)averaged quantity denoted by a tilde over the variable is 
define as 


20] tothe af 
—— dt = -— 21- 
sf 554-5 (21-302) 


where, as before, a bar over a variable indicates a Reynolds-averaged quantity. It 
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is important to note that Favre-averaging is a mathematical concept introduced to 
simplify the equations, and it is not a physical phenomenon. 
Relation (21-302) can be written as 


pf =of = (+e) (F+f) =P f+ ef (21-303) 
Any quantity f can be decomposed as 
f=f+f" (21-304) 


where f and f" are the Favre mean and Favre fluctuating part of the turbulent 
motion. Relation (21-304) can be multiplied by p and time-averaged to yield 


pf=pf+pf" (21-305) 
Recall that from relation (21-303) 
pf = pf (21-306) 


and therefore it follows that 
pf =0 (21-307) 


Now, from definition (21-304), one has 


f'=f-f (21-308) 
and from (21-303) E 
f=f+ ar (21-309) 


The combination of relations (21-308) and (21-309) yields 


"Ro Pe api Pf ,_ PF M 
F-f-f-f-f F f 5 (21-310) 


where relation (21-19) is also used. Subsequently, the time-averaging of (21-310) 
provides 


m.p Lf 
F=f 7 


But since, according to (21-21), ff = 0, then 





Fi WEL #0 (21-311) 


Turbulent Flow and Turbulence Model: 10 


Now, the rules of Favre-averaging are summarized as follows 


f=f (21-312) 
pf" = 0 (21-313) 
eiu. EHE f 

f ; (21-314) 
pig = pÍóvpPg (21-315) 


With the definitions and rules of Favre-averaging established, consider the time- 
or Reynolds-averaged continuity equation (for a two-dimensional flow) as given by 
Equation (21-23), that is, 

Op O yo — ð s Ss 

aàr* a; U pü gu) s. (p -- pv) =0 (21-316) 
Using relation (21-303), the Favre-averaged continuity equation now can be written 
as 

ô ð ð 

= + x (pt) + — (07) =0 21-317 
DE t a; PË + a, 09) ( ) 
Observe that the Favre-averaged continuity equation retains the same form as the 
Reynolds-averaged continuity equation, and, in fact, the same form as the original 
continuity equation. This is precisely the reason for the definition of mass-averaged 
quantities, as given by (21-303). That is, the mass-averaged quantities are defined 
in such a way as to write the resulting equation in a similar form as the original 
equations, while retaining the effect of density fluctuations. Furthermore, this is 
the reason why it was stated previously that the concept of Favre-averaging is 
mathematical. 

Now consider the z-component of the momentum equation for a two-dimensional 
flow given by 

ð à cm OT. OT xy 
z (pu) + 5z (pu?) + 5 S Mays i ne Dr + du (21-318) 

For simplicity, each term will be considered separately. The procedure, as illus- 
trated previously, begins with time-averaging, and subsequently the Favre-averaging 
concept is applied. Now the first term is written as 


209-2 [GFE] = 2 es 7v) - 2 G4 


Similarly, the remaining terms are modified and given by the Lun 


Ó —. ð .2,—z 
a; Pw) = g PP pu") 
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PEE E Mr 
ay 97 = ay Piet ew v") 
noo (49829 
= = F\3 a2 30y 
and 


E .[0ü , dv 
Try — (= + a) (21-319) 


Therefore, the z-component of the momentum equation, given by (21-318), is writ- 
ten in a mass-averaged form as 


2 pa? gr? (pai) = -2,2 (a DEES Cy pura) (21-320) 


Note that the density in pu” and pu" v" is the instantaneous density, and that 
these terms represent Reynolds stresses. 
Similarly, the y-component of the momentum equation is written as 


P092 (pii) z 09) = -2+2 (Fay — pu v5 = y; cp?) (21-321) 


The energy equation in terms of the total energy per unit mass expressed in Carte- 
sian coordinates for two-dimensional problems is given by 


a a a 
ay 099 + az; (ouer + pu) + y (ove, t pv) = 


x (U Trz + v Tz, — qz) t» 2 a Tay + U Tyy — Qy) (21-322) 
Again, considering the modification of each term separately, the first term is 
& 097 5 [GF MED] = AE 5 (pe) 
The term pé, can be rewritten as 
pi =pG=paTt+= CEN) (21-323) 
where the following relation is used 


1 1 
& — et 2V - eT S (uv) 


2 
Observe that, in the equation above, the potential energy has been dropped. 
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Now the instantaneous terms in Equation (21-323) are replaced by their mean 
and fluctuating values, thus, 


— Hun E e 
P&T + 5 p(u* +o?) e (P+ e)(T- T) 5 pur+ pv? 


MEE, NE Rer 
7 e(pT * eT) +> put P p 


Using relations (21-309) and (21-315), one can now write 


2a deus 1 lc 
LENS oT 24d) Lou" .LlLn 4 l5 
p& = «(p )*t3pü + spun + opi pv 
or 1 1 
Pë = e pT + Pl +0) + 5 (pu? + pv?) (21-324) 


The last term in relation (21-324) is defined as the turbulence kinetic energy (per 

unit mass) as follows 
1 
pk Lp Ev 
or 
1p v5 
2 p 

It is important to note that the turbulence kinetic energy defined in the Favre- 
averaged formulation is different from the turbulence kinetic energy defined in the 
Reynolds-averaged formulation given by (21-75) as 


k (21-325) 


ks ju +0) (21-326) 


which is essentially based on incompressible flow considerations. Now relation (21- 
292) can be written as 
Pè =e pT + =p (a+ 6) + pk (21-327) 


The second term of the energy equation (21-322) is considered next, which is rewrit- 
ten as 


1 
pue puc (pes pu [p T z p+) e p] (21-328) 
Using the perfect gas equation of state, one may write 


pe, T-- p — pe, T t pRT = pT(c, 4 R) - pe,T 
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Thus, (21-328) is written as 
1 1 
pco T + g^ v) u —cpuT-* ; pu(u +v) 


which can be time-averaged to provide 


pul = püT + pu T" 
according to (21-315), and 


1————— l; lL 
zprvu + v?) = ze% t qeur 


1 — — 
5 pü(ü^ 4-0) +a EIS n v?) 4 üpu^ 


zs | 
3j Wy 3 rl 
+ op uv" + 2 pu" + pul!) 


1 — ee 
x pi (a? 4-0?) - pkitipu” + pu" 


1 puu) 
+ 2 pu” + pully’ 


Finally, the second term is written as 


———— B 1 Me SE 
(patp)u = pu ot +5@ +0) ek] es p" 


—— , 1 
+ iipu” +o pulv? + 3 (pu + pul”) 


Recalling the relations 
=p 1 afzI2 ~2 1 T P\ a2 
e, pT + z p? +0) +5 (p"  pu'v^) = pe 


and 
where 


Then (21-331) is written as 


(pet pu = (Pë +p) ü+ ceppu T" + üpu^ -- üpu'v" 


+5 (pu + pur) 


(21-329) 


(21-330) 


(21-331) 
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Similarly, the third term is modified to yield 


(pe pv = (Dë +p) e pv" T" + ü pu" +i pv? 


1 
= m " 
+5 puu! + pu) 


The terms on the right-hand side of the energy equation include shear stress terms 


such as U Tz,. These terms, for example, are written as 
Utz: = (ü + u”) Tax = Ü Fer + U" Tre 


Finally, the Favre-averaged energy equation is written as 


ô à ee lI 
a; (Pe) + az [Gà 5) Ü+ c, pu" T" + ü pu” + üpu'v" 
1 B L puihi NN Say. se r 
+5 pu + pur) | e x. [Gà B) 5+ e p T? + ü pu 
5 1 1 IP gyi "3 à ? T. ulr 
TÜpv +5 (pu^ v + pv ) ns U Tre + U" Tez 


Big Mba On lee see ee a eee m 
+Ü Tay + U" Try — Ve lta [tay + Wray eS + Phy d, | 
and rearranged as 
Qu. Ô som aa LT 
pe Pa) + a; Le + Pa) + s [Ge +9) 5) = 
ð y RS Mal me = es 
óz [-àpu —UVpu"v + Ü Far + Ü Tay] 
C a8 zz ane 
+ oy |-üpu —ÜUpv uL ELM 
à —J— 8 Li 
+z -opu — ae] + ac [-e pv"T" — 4, | 


à 1 — fy oe 
+ n l-3 pu” + p uv" ) + Ten + vs 


ð 1 TAE 
TE l-5 pu v" + pv! + Wig t | 
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The terms involving triple products such as p u^, p uv" , etc. are typically smaller 
than other terms and are usually dropped. Therefore, the energy equation is reduced 


* 54 )] 


(21-332) 


to 
2 pa) + 2 (Garp a+ È (08 Dil - 
x [a (pu + tex) + 9 (-puv + Fey ) | 
+e fü (PU + Fey) + (-pv™ 
+2 Lepr -a]+ > [opra] 
where 
pop. 4.107 
& = -k gz pr ðr 
and 


Now the system of equations, composed of the Favre-averaged Navier-Stokes 


equations, are written in a vector form as 


Q , OE | OF _ ab 


where 


(21-334) E= 


üt Oe Oy ^ Os Oy 


(21-333) 
 pü 
UN (21-335) 
pid 
(p& + p)ü 
(21-336) 
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0 
Tox cu pu” 
E, = E Lut (21-337) 
Fry pulv 
di (fs "Y pu^) 4 (Fey -7 uv) — d — Cp pul T" 
0 
Fey — pulv" 
E o- 21-338 
v Ty _ p v^ ( ) 


ü (Fey — puo) +6 (Fy — pv?) — a, — c pv" T? 

Similar to the Reynolds-averaged equations, the Favre-fluctuating terms are re- 
lated to the Favre-averaged mean quantitites through the introduction of turbulent 
viscosity. The procedure and definition of terms are provided in the following sec- 
tion. 


21.8.1 Turbulent Viscosity 


Following the Boussinesq approximation introduced in Section 21.3.1.1, the 
turbulent shear stresses are related to the rate of strain in a similar way as the 
laminar shear stresses. Thus, a turbulent (or eddy) viscosity is introduced in the 
expressions for turbulent shear stresses as follows. 


gu m (5 " 5) (21-339) 
-PF = m (25-30-07) ~ 2p (21-340) 
Mor eee (2-iv 7) -25k (21-341) 


The stress terms defined above can be expressed in a tensor notation given by 
—+ FF Oü; 84; 20ü, 2 

"gb — nd hou! ead ee ee E a se » 

ds (5: Or; 30, is) 3 pkóy 

Relation (21-342) can be written in the following form as well 

1 die 

3 Ox; 


(21-342) 


oe " 2 
— puju] = 2j (3 = is) — g bk by (21-343) 
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where 


[e : (= à as) (21-344) 


Turbulent heat flux terms are also written in a form such as to resemble the laminar 
heat flux term. For this purpose, a turbulent thermal conductivity is introduced 
such that 


m oT 

c, PUNT” = —k (21-345) 
and $ 
PEU oT 

c p v" T" = —k, By (21-346) 


It is a common practice to write the turbulent thermal conductivity in terms of 
turbulent viscosity by the introduction of a turbulent Prandtl number defined as 


This is the same approach taken in the Reynolds-averaged formulation. Typically 
the value of the turbulent Prandtl number is selected as a constant, and for air is 
specified as either 0.9 for wall bounded flows or 0.5 for free shear layers. 

Now the stress and heat flux terms in (21-337) and (21-338) are redefined as 
follows. 


Tac Tyr — pu? = (+ m) |2 J% -— 20. ğ |- 2 5k (21-347) 
ôr 3 3 
— Oi | O0 
= E NEUES T — (pn balers E E 
Try Tzy puv (i + He) (5 + x) (21 348) 
Ty = Ty pu = (ji + m) ee wy SV 2 -— pk (21-349) 
, "rc oT T aT 
se = rae p = (kka) dz (5 $T Pr) E (21-350) 
, —— gh oT oT 
-a = -p-p GR) =o (f+ pe) a CD 


Finally, the FANS equations are written as 


0Q OE OF _ OE, F, 
Ot tr z d Oy ðr - *X 








(21-352) 
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where 


dD DI 
Re 
DI DW 
& gm 
+ 
"3r 


Q- (21-353) E- (21-354) 


DI 
Qu 
DI 
Rr 
e 


DI 
Qi 
pom 
"OI 
e 
+ 
"3 
— 
g 


(21-355) 


E, = (21-356) 


FA = 7 (21-357) 
Tyy 
ü5. + Fy — d, 


Observe that the flux vector formulation given by (21-352) and the flux vectors 
E, F, E,, and F, given by relations (21-353) through (21-357) are in a similar form 
as the corresponding terms in the Reynolds-averaged formulation, and, indeed, in 
the same form as the corresponding laminar formulation given by (11-192) and flux 
vectors (11-193a) through (11-193c), (11-193e), and (11-193f). 


21.9 Concluding Remarks 


Turbulence is a very complex physical phenomenon. Direct simulation of the 
governing equations with small scale turbulence for a general configuration at high 
Reynolds numbers is not currently possible. The available approach for practical 
applications, at the present, uses various types of turbulence models which relate 
the average effect of turbulence to mean flow. A tremendous amount of work is 
still required in turbulence modeling. Currently, a universally accepted turbulence 
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model which can be used in diverse flow regimes does not exist. Some models 
may perform reasonably well in certain flow regimes but result in inaccurate pre- 
dictions under different situations. Therefore, turbulence models must be carefully 
investigated for their range of applicability before being utilized. 

An attempt has been made in this chapter to introduce some fundamental con- 
cepts of turbulent flow, turbulence models, and limited number of solution schemes. 
It is hoped that this attempt will be beneficial as a first step toward the inclusion 
of turbulence into the governing equations of fluid motion. An in-depth review of 
turbulence may be found in various texts such as [21-25] through [21-27]. 
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22.1 Introductory Remarks 


In order to numerically solve a partial differential equation, the differentials 
appearing in the equation must be approximated by algebraic expressions. Taylor 
series expansions are typically used to develop these expressions as illustrated in 
Chapter 2. Several formulations expressed as forward, backward, or central dif- 
ference approximations with various orders of accuracy can be easily developed. 
Typically, the expressions for the partial derivatives are developed and expressed 
solely based on the dependent variable. For example, a second-order and a fourth- 
order central difference expression for the first derivative are given respectively by 

Of ^», fa-fa 2 
> =f,= Az + O(Az) (22-1) 
and 


Of) p fea ~ 8h + Bf = fus 
ðs C' 12Az 


The grid points involved in the computations are shown in Figure 22.1. 

Essentially what is being accomplished in the estimation of the derivative f at 
point i is a curve fit through several points of data and subsequently determination of 
the derivative. In fact, the finite difference approximations can be directly obtained 
by a curve fit of polynomials as shown in Chapter 2. 

Next logical approach to increase the accuracy of the estimates of the deriva- 
tives, in particular for problems involving shorter length scales or equivalently high 
frequencies is to include the influence of the derivatives of neighboring points in 
the calculations. However, keep in mind that these "neighboring" derivatives are 
unknowns, and therefore the developed expression will involve more than one un- 
known. Nonetheless, when this equation is applied to all the points within the 
domain sufficient equations are produced which can be solved simultaneously for 
the unknowns. This approach is analogous to the solution of a pde by an implicit 


+ O(Az)4 (22-2) 
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scheme (more than one unknown per equation), compared to an explicit scheme 
(one unknown per equation). The resulting approximation is referred to as Com- 
pact Finite Difference (CTFD) formulation. They have been also referred to as 
Hermitian formulations. 'The derivation of the CTFD formulation also involves 
Taylor series expansions. À simple example is used here to illustrate the procedure, 
and subsequently a general formulation is provided. 








Figure 22-1. Grid points appearing in formulations (22-1) and (22-2). 





22.2 Compact Finite Difference Formulations for the 
First-Order Derivatives 


Consider a three-point Hermitian formula involving the first derivative given 
by 


m=t1 


Hi = Y (am fim + bmfitm) = 0 (22-3) 
m--1 
or 
fiai aofi + aifi- + bifa + bof; + b-af1 =0 (22-4) 


Substition of Taylor series expansions into Equation (22-4) yields 


si 1 2#, 1 ae, 1l agw, l 5 
ai [fit Arh + (Az) fi (A2 KE gA fi +A E 
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zi; (89 ff + (Aa ] + aofi + a-i [f - Aah + (88 ££ - KA2? + 
1 iv 1 i 1 = t " 
gi (A9 ff - 15 (A ft + aon (A2 ff O(Az)] b [fi nt 


1 


1 27" 1 3 riv 1 4 fu 
SAP + ZAIE 2 (Aa + i 


(Az) f?! +O (Az) x 


^s d 1 2 » 1 3 riv a. Arpv.— 
b, f Az f, + 5 (Az) fi g (As) i + 54 (42) fi 


EN 


gg A9 fi + O(Az)*] =0 (22-5) 


or 


(a1 + ao + a1) fit (mi — ai) Az +b + bo + dal K+ | (hm ( +01) + 


Az(bs — 0-2] ft | s (Az) (a — a=) + F(A)? + ba] A + 


| g; (A9 +aı)+ s (A20 - b..)| ff + | i; (40a —aı)+ 
5, (A9*( t m R + | z;; A2 + a) + 
sy tabi -b| - o (22-6) 


Now if the relation given by Equation (22-6) is exact, then the coefficients in (22- 
6) should be zero. In reality, however, only a few of the coefficients are set to zero, 
and the remaining higher order terms are truncated and will form the truncation 
error (TE). To obtain a third-order scheme, set the coefficients of f;, Åi fis and f, 
equal to zero. 


01 * a9 +a — 0 (22-7) 
Az(ai = a) +b + bo + b_1 =0 (22-8) 
a (Aa) (a. - a1) + (Az) — 6-1) =0 (22-9) 
s (Az —aa)-c NA +b_1) =0 (22-10) 


The coefficients are a), ao, a1, and bg are solved in terms of bi, and b; to yield 


1 
Q1 = zAz C 9b - b-1) (22-11) 
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2 
a = Az T bi) (22-12) 
1 
ay; = Dag (t 9b-1) (22-13) 
bj = b o3) (22-14) 


The truncation error from (22-6) is 


TE - [Ara E s (Az) (b 3 2] fü EOG any 


1 4 v 1 6 1 5 vi 
aAa) b tb] A+ [s (hm (ie aa) + (An) — b4)] A+ 
or with the substitution of (22-11) through (22-14) 
TE = Z (Ar) — ba) + sc (Az) (b t b) fi? zac (Am) (bi — baa) £i" 
12 Ve 60 ~*~ 180 virt 
Te (22-15) 
The formulation (22-4) can be written as 
ay do 0 ' bo pt bap 
EAA + bi fit p + fit b, f+ b. fa (22-16) 


where the equation is divided by bj. Now, the expressions given by (22-11) through 
(22-14) are redefined as 


ay, 1 

ZE ee) eee (2 Dos 22-1 
5 2A; « oO (22-17) 
a 2 

— = (1- 22-1 
b, DA (22-18) 
a_; 1 

LOE tes ME OE 22-1 
m ZA; * 5e) (22-19) 
Pee, ae (22-20) 
bi 

where a = bb. 


Substitutions of (22-17) through (22-20) into (22-16) yield an equation for the 
derivatives f; ,, f,, and f,,;, as follows 
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a f, .,4-2(1-- o) f Ju = -gA; (1450) f- =l- o) fit a sig f TE 
(22-21) 
where 


TE- ij (As) - a) fi? + aj (A2*Q + a)fi + —- (Az) (1 — a) fi +.. (22-22) 


TS 
Note that for œ = 1 the leading error coefficient vanishes and the scheme becomes 
fourth-order, whereas, for other selections of a’s the scheme is third-order. For 
a = 1 the formulation is 


Ratan + fin = m Ac fi) (22-23) 


Now a general five-point formulation for the approximation of the first derivative 
can be written as 


4 , " n ' af fi-1 fi-2 
Biata hiit fi tafi tB fia = ery + Vet fen Be 


fas — fi-s 
6Az 
(22-24) 
The relations between the coefficients in Equation (22-24) are established by 
matching of the coefficients obtained by the substitution of the Taylor series expan- 
sions. The procedure is illustrated for the second-, fourth-, and sixth-order schemes. 
By retaining higher order terms in the Taylor series, the relations for the eighth- 
and tenth-order schemes can be established. 
The Taylor series expansion of the terms on the left-hand side of (22-24) is 


4 " " 4(Ax)? m 8(Az) 

















fia = f-üAmh ok - At + IY fe (22-25) 
fia = -Asfi + , G2 K- Cay A + en f (22-26) 
fin = nume Ln i Mf (22-27) 
foe = f+ 2dafl + SOR pry sary AAT fwg [ov DU f? (22-28) 


The left-hand side of (22-24) is 


B fist ofi fib odas B fuam Bra fas) t o(f ei fen) + f. (22-29) 
Substitution of expansions (22-25) through (22-28) into (22-29) yields 
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32(Az)" , 


plog e Aon + OE ne] + alag + nn? fe mx 


4! 








J+ f= 





(2a + 28 + 1) f; + (Az)'(a + 48) f. + (Gay (a+168)f,° (22-30) 


The required Taylor series expansions for the terms on the right-hand side of 
(22-24) are 


























fies = fic sf; + WARM pr PERN pe q SAY po _ MX e (22-31) 
fia = fi-2Aaf + Kas gm e, was) a - FAA) e (22-32) 
fin = h angia OED g (OOP pr | (Ba py (BY py TE 
dave d as; ga oe P" Gm (Az pus 4 E Me (22-34) 
fas = fit 2Azf, + T" —— + Rs pe 1A ya 80 e (22-35) 
faa = fit 3AF; + xa fe LOR ge, AT pe 2308. e (22-36) 


Expansions (22-31) through (22-36) are substituted into the right-hand side of 
(22-24) to provide 


zz (Fest — fii) + garlar- fea) + gclos- fis) = 














A A 3 
is | 2o ex s fex s si] a 4Ax [Aa fi rem f 
A 5 3 5 
2] "£x feaz f+ On A pias) ; 


Cay" 








= (a+b+c)f; + (a+ 454-96) f +4 Ary (a+ 16b + 81c) f? (22-37) 
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Now from expressions (22-30) and (22-37) one obtains 


(20 + 28 +1) f; + (Az) (a+ 46)f," + `- : 
(Az)? 


3! 





(a + 168) f? = 


(Az) 


(a+ 4b--9c)f, + : 


(a--b--c)f, + 








(a + 16b + 81c) f; (22-38) 


Setting the coefficients of the differentials equal, the following constraints be- 
tween the coefficients are established 


Second-order: 1+2a+28=a+b+ce (22-39) 
! 
Fourth-order: 25 (a+ 228) = a+ 22b +3% (22-40) 
! 
Sixth-order: 25a + 218) =a + 21b + 3c (22-41) 


The constraints for higher-order terms are : 


[| 
Eighth-order: 27(o +288) = a+ 255 + 35c (22-42) 


Tenth-order: 25 (o 4-288) = a+ 25b 4- 3*c (22-43) 


To obtain the error term of a particular scheme, the truncated terms are con- 
sidered. For example, for the fourth-order scheme, the truncated terms are 


TE- [g(a + 2b + 3%) — Zla + 248) (Az) f; (22-44) 
For the values of 8 = 0, a = ¿(a + 2), b = &(4o — 1), and c=0 


TE = gj Go - D(Az)*f? (22-45) 


Depending on how aand f are specified, either a tridiagonal system or a pen- 
tadiagonal system will be generated from Equation (22-24). For the special case of 
B = 0, the resulting system of equation will be tridiagonal. Several special cases 
are presented in this section. 

For B = 0 and c = 0, a one-parameter (in ternis of o) family of fourth-order 
tridiagonal systems can be produced where a — 2(a--2) and b = (4o — 1). As seen 
previously, the truncation error of the scheme is T'E — & (3o -= D(Az)!f; 


Observe that for a = 1, then b = 0, and a = 3, and one has 
ly ro le _3 fin-fina 
gfi- + fe + qin =3 Skane (22-46) 
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which is identified to Equation (22-23). Furthermore, note that the selection of 
a = 1/4, b = 0, and a = 3/2, satisfies the requirement given by relation (22-40). 
The one-parameter family of schemes are summarized in Table 22.1 











Order of 


a 
accuracy 


je9 | 42-9 | 


1 (568a— 
1(12—7a)| 30 5 





Q 
. Max LHS 
22-47 £(3a — 1) (Az)*f? 
2 (—8o + 3) (Az)6 f?" 
M4 (20 — 1) (Az) f? 


Table 22.1. Selected one-parameter family of compact finite difference schemes 
for the first-order derivative. 





The following observations can be made with regard to the schemes which are 
given in Table 22.1: 
1. Scheme I 
(a) For a = 0 the values of a and b are $ and -}, respectively, and relation 
(22-24) provides 
pa fi-2 — 8fi-1 + 8fit1 — fire (22-50) 


12Az 
which is the fourth-order central difference formulation. 


(b) For a = land the resulting values of a = 5 and b = 3, the scheme 
becomes formally sixth-order due to vanishing of the leading term in the 
truncation error. 

2. Scheme II 

(a) Fora =i, then a = &, b=, and c= 0, and the result is identical to 
formulation of Ib. Therefore, the formulation of Ib can be considered as 
a member of II. 

(b) For a = 3, then a = 8, b = }, and c = —g;, and the scheme is eighth- 
order. 

3. Scheme III. 

(a) For a = $, p =0, a = , b=}, c= 

(b) For o = 4, 2 = $, a = Ë, b = M, and c = y, the scheme becomes 
tenth-order accurate. 


—g the scheme is identical to IIb. 
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22.2.1 One-Sided Approximations 


Several schemes identified for the first derivative based on the central difference 
approximation of (22-24) will encounter difficulties for nonperiodic boundary condi- 
tions. Typically, one-sided approximations are used at the boundaries. Considering 
points ? = 1 and i = IM to be the boundary points, the following general expressions 
for the first derivative may be written. 


fi + of; = aah tbh + cfa + dfa) (22-51) 
and 


t + 1 
fim + afim = Az (fiM +bfimm + Cfimme + dfimm3) (22-52) 


The relation between the coefficients are established in the same manner as pre- 
sented previously and are given by the following. 


Second-order: a= — 568 +a + 2d) 
b = 243d 
c= - la. o. 6d) 
2 
and the truncation error is 
1 " 
TE = a — a ~ 6d)(Az) f, (22-53) 
: 1 
Third-order: a= — gt + 2a) 
E ;(6 ud 


1 
1 
d = sa = a) 
with the truncation error of 


TE= BC — 3) (Az)? i (22-54) 


Fourth-order: a 


ll 
w 


a= —— 
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and 


TE = (Ar) fe (22-55) 


22.3 Compact Finite Difference Formulations for the 
Second-Order Derivatives 


A similar expression to that of (22-24) can be written for the second derivative as 
follows 


Bap AG aa Bh ee pe M den pfe — 2A + fia 


(Az? (Az? 
fas — 2fi4- fi-a 
dec A (Ane (22-56) 


The relations among the coefficients are established as 


Second-order: 2a+284+1 = a+b+c (22-57) 
! 
Fourth-order: je +278) = a4 23b 3c (22-58) 
! 
Sixth-order: qe +278) = a+%b+ 3'c (22-59) 
! 
Eighth-order: ge 4288) = a+ 2%b+ 39c (22-60) 
10! M 
Tenth-order: ge 4-298) = a+ b+ 2%c (22-61) 


Several formulations for the second derivative are provided in Table 22.2. 
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č Order of 
accuracy 
4th 


$(10a - 1) | 4h | 
1(6-9a |1(-3--24e| i(2- 11a "m 
—128) —69) +124) 
zu (38a | 5 (696 2 (245a | lla ae 
-1191a) —294) —344) 


—8(9 — 38a + 2148) (Azx)® f$ 5 
zam (8990: — 334) (Az)? f* 5 


Table 22.2. Selected compact finite difference schemes for the second order 
derivative. 





Observe that for the scheme given by Equation (22-62), when a = 0, the values 
of a and b are determined to be $ and —1i respectively. Subsequently Equation 
(22-56) provides 


"o Ihat fir- 30fi + lfi- fir 


Í IXAzyi (22-65) 


a=2,a=#, andb- $, the scheme is formally sixth-order accurate. 


22.4 Error Analysis 


To investigate the resolution characteristics of finite difference approximations 
and to perform error analysis, the Fourier analysis approach is undertaken. This 
procedure is relatively simple and provides an effective method to quantify the 
resolution issue of the finite difference scheme as well as error analysis. 

To further simplify the Fourier analysis of error, a function which is periodic 
over the domain is assumed. Now, the function is expressed with a complex Fourier 
series by 


k=IM/2 
fz)- X Rev (22-66) 
k= -1M/2 


where I = /—1, L is the length of the domain, z is the independent variable, and 


F; is the kth component of complex Fourier coefficient. Furthermore, Az = iz is 
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the spatial grid size, s = à and K = Zz = ot are defined as the scaled wave 
number. Now, 
k=IM/2 
fa y. Ent (22-67) 
k=-IM/2 
and the first derivative of (22-67) is 
: k=IM/2 
f()- Y FIK e” (22-68) 
k=-IM/2 
Or f k=IM/2 ; 
f()e X Re" (22-69) 
k--IM[2 
where 
F,-IKF, (22-70) 
Now, consider a finite difference approximation of the first derivative of f as 
: IM/2 ; 
fus)- Y; Eue” (22-71) 
k- -IM[2 
or nis , 
fals) = Y, RIKe*' (22-72) 
k=-IM/2 
where , "i 
FQ-IKEF (22-73) 
and K' = K'(k) is the modified wave number. One obtains the exact value if 
Fy, = Fy (22-74) 
or 
k'=K (22-75) 


which would be a straight line in a wave number, modified wave number coordinate. 
Now, define the error as 
_K-K i 





ER = K (22-76) 
from which one may write 
K' = K(1— ER) (22-77) 
From the definition given by (22-76), the following may be written 
‘ EY y IKs _ !IKs 
PRE f (s) — frals) _ DELK e EFIK e (22-78) 


f'(s) E FLIK elK* 
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Now, consider an example where a sine wave is used. In this case, the Fourier 
series is reduced to one term only, therefore, 


f(z) = sin(Ks) = sin(K ©) (22-79) 


and 
D K gz K 
f (=)= Az (K AD me cos( K s) (22-80) 


APPROXIMATION OF FIRST DERIVATIVE 
—~ Exact 

—— Compact finite difference (4th order, alpha=5/14) 
—¥— Compact finite difference (6th order, alpha=1/3) 

—©— Compact finite difference (4th order, alpha=1/4) 

~—A— Central finite difference (4th order) 


E 
bæ 
& 
E 
3 
G 
m 
> 
Ld 
= 
v 
£ 
hz 
5 
© 
z 


1.50 
Wave number K 





Figure 22-2. Comparison of several schemes for the approximation of the first 
derivative of (22-79). 


A comparison of the wave number K versus the modified wave number K' is 
shown in Figure 22-2 for several schemes. It is evident that the compact finite 
difference follows the exact differentiation over a wider range of wave numbers com- 
pared to the fourth-order central difference approximation. The approximate value 
of maximum wave number for well resolved waves (k.) is provided in Table 22.3. 
À measure of resolving efficiency of a scheme may be defined by e = Be For ex- 
ample, as seen from Table 22.3, the resolving efficiency of the fourth-order scheme 
with a = 5/14 is higher than the resolving efficiency of the sixth-order scheme with 
a= 3 
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Scheme 


Central finite difference (fourth-order) 


Compact finite difference (fourth-order, a = 1) 


Compact finite difference (sixth-order, a = i) 


Compact finite difference (fourth-order, a = £) 


Table 22.3. Comparison of the resolving efficiency of several 
Schemes. 





22.5 Application to a Hyperbolic Equation 


A simple linear hyperbolic equation is used to obtain the solution by several 
schemes and compare the solutions to each other and the analytical solution. The 
numerical schemes considered are (1) a first-order upwind scheme, (2) a modified 
four stage Runge-Kutta scheme where the convective terms are approximated by 
a fourth-order central difference expression, and (3) a modified four stage Runge- 
Kutta scheme where the convective terms are approximately by a sixth-order central 
compact finite difference scheme. 

Consider the first-order wave equation 
oe = -at a>0 (22-81) 
where a is specified to be 40r m/s, and the domain of solution is 0 < z < 10r m. 
An initial wave is specified as 





u(z,0) = 0.0 0zzr«2m- 
u(z,0)- 1.0 — cosx 2r X z « 4v (22-82) 
u(z,0) = 0.0 4r € z « lOr 
The spatial grid is 
Ar= 107 
IM —1 


where IM = 51. 
Several temporal steps (or corresponding Courant numbers) are specified as 
follow. 
1. (a) At = 0.005 sec, c=1.0 
(b) At = 0.0025 sec, c=0.5 
(c) At = 0.00125 sec, c= 0.25 
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The explicit first-order upwind scheme applied to Equation (22-81) yields the 
following FDE. 


ut = uf efu- ut.) (22-83) 
The modified Runge-Kutta scheme is given by 
uf) = ys (22-84) 
D . Q»& At dua) f 

ui za i a 4 az (22 85) 

At ðu 
uf) = u- ayla P (22-86) 

n At ðu 
uf = u “ea 9 (22-87) 
ur = yt aati (22-88) 


Two options are considered in the approximation of the first-order derivative 
on the right hand side of Equations (22-85) through (22-88). In the first option, a 
fourth-order central difference approximation given by 


Ou _ uia —8u,. + 8ui1 — uiia 4 
A UTAH + Q(Az) (22-89) 


is used for 3 < i < IM — 2, where a second-order scheme is implemented at i = 2 
and i = IM -— 1. 

For the second option, the derivatives are evaluated by the central compact 
formulation given by Equation (22-48) with œ = 1/3, a= 14/9, and b=1/9, 
which is formally six-order accurate. 

Before the inspection of numerical solutions, a review of several observations 
with regard to the numerical schemes to be used is in order. 


l. Recall that the error term for the explicit first-order upwind scheme is given 
by Equation (4-77), repeated here for convenience. 





ER - a - oot - LAT og -3c4- yes + 
O[(Az)’, (Az)(At), (Az)? (At), (Azy]] (22-90) 


Note that the dominant term involves a second derivative 02u / Oz? , and there- 
fore the error produced in the solution is dissipative. That is, as the Courant 
number is decreased from its maximum value of one, the amplitude of the 
wave is decreased and the wave is dissipated to the neighboring points. The 
solution at Courant number of 1 is exact as seen from Equation (22-90), where 
the error becomes zero. 
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2. The dominant error in schemes which utilize central difference approximation 
for the convective term is dispersive and results in oscillations in the solution 
in the vicinity of large gradients. These observations are clearly evident in 
Figures 22-3 and 22-4. The solutions by the first-order upwind scheme for 
the three different Courant numbers (or time steps) are shown in Figure 22-3. 
The solution for Courant number of 1 is exact as seen in Figure 22-3a. The 
dissipation error and its effect on the solution is shown in Figures 22-3b and 22- 
3c. The solutions by the modified Runge-Kutta scheme with central difference 
approximation of the convective term are given in Figures 22-4a through 22- 
4c. Several schemes are available to reduce the oscillations, for example, 
addition of a fourth-order damping term or incorporation of a TVD model into 
the formulation as shown in Chapter 6. The solutions by the Runge-Kutta 
scheme with central compact formulations are illustrated in Figure 22-5. The 
solutions retain the original wave shape for all Courant numbers specified in 
the problem. A comparison of numerical solutions with each other and the 
analytical solution are shown in Figure 22-6. 


The solution by the compact formulation is well behaved for all the three spec- 
ified Courant numbers. Of course, this increase in accuracy is accompanied by an 
increase in computation time. However, it should be stated that in some applica- 
tions it will be necessary to implement the compact formulation due to its resolution 
efficiency of small scales, which may be present in these applications. Example of 
such applications includes problems in computational acoustics and direct numerical 
simulation. 
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Figure 22-3. Solutions of the wave equation by the first-order upwind scheme. 
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Figure 22-4. Solutions by the fourth-order Runge-Kutta scheme with central 
finite difference approximation of the convective term. 
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Figure 22-5. Solutions by the fourth-order Runge-Kutta scheme with central 
compact finite difference approximation of the convective term. 
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Figure 22-6. Comparison of the solution of the first-order wave equation by 
several schemes. 
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22.6 Problems 


22.1 Consider a wave which is propagating within a domain with the following 
initial distribution and shown in Figure P22.1. 


u(z,0) =0 0cz 2m 
u(z,0) — 0.5 {sin [az H 77)/2 ü 1} Qn < z < 81/3 
u(z, 0) = sin [ss - 2)/2] 81/3 < £ < 107/3 
u(a,0)=05{sin|(@e~11n)/2]-1} 107/352 <a 
u(z,0) = 0 4r < z € 10r 





-—Ó—— PER eA 
15.0 i 25.0 
x 


Figure P22.1. The initial wave distribution at t = 0.0. 





The governing wave equation is given by 


Lr ce 
8t 


where a is specified as 407 m/s, and the domain of solution is 0 € x < 107. Select 
a spatial step of Az = (107/IM — 1)m where IM = 51, and temporal steps of 
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(a) At = 0.005 sec, c=1.0 
(b) At = 0.0025 sec, c=0.5 
(c) At = 0.00125 sec, c= 0.25 
Obtain the solution up to t = 0.1 sec by the following schemes: 
(I) Explicit first upwind scheme 
(II) Modified fourth-order Runge-Kutta scheme 


(III) Modified fourth-order Runge-Kutta scheme with central compact finite differ- 
ence approximations 


Plot the solutions at intervals of 0.02 sec up to 0.1 sec for all three schemes and for 
the three Courant numbers. 


22.2 Repeat Problem 22.1 for the following initial velocity distribution given by 


u(z,0) = 2exp - Eal 


and shown in Figure P22.2. 





Figure P22.2 The initial wave distribution for Problem 22.2. 
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Large Eddy Simulation and 
Direct Numerical Simulation 





23.1 Introductory Remarks 


Computation of turbulent flows has always been challenging and remains so at 
the present and in the foreseeable future. However, since most applications involv- 
ing fluids include turbulence, a push toward the solution of turbulent flow has been 
a continuous endeavour. As discussed previously, essentially three avenues are avail- 
able. First isthe solution by the Reynolds Averaged Navier-Stokes (RANS) equation 
complemented by a turbulence model. Several turbulence models are available for 
this purpose. These models are continuously evolving (modified and improved) in 
order to increase their accuracy and range of applications. However, since all scales 
of turbulence must be modeled in the RANS equations, a model which is capable 
of predicting turbulence over a wide range of flow conditions and geometries has 
been illusive. Nonetheless, the RANS approach for computation of turbulent flows 
is very practical, and the resulting solution is relatively accurate if an appropriate 
turbulence model is used. 

The second approach is the Large Eddy Simulation (LES), and the third ap- 
proach is the Direct Numerical Simulation (DNS). Application of DNS for flowfield 
analyses and design purposes at the present is not practical due to limitations on 
available computers of today and due to the enormous amount of time required. 
Large Eddy Simulation requires less grid points and computation time, typically 
five to ten percent of CPU time of DNS computations. However, LES computation 
still requires substantially more CPU time than similar RANS computation. LES 
and DNS solutions are, however, very valuable as they can be used to improve tur- 
bulence models for RANS and to gain better understanding of physics of turbulence. 
The objective of this chapter is to introduce LES and DNS. 
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23.2 Large Eddy Simulation 


Large Eddy Simulation provides a compromise between DNS, where all scales 
of turbulence are computed directly from the Navier-Stokes equations, and RANS 
equations, where all scales of turbulence must be modeled. In Large Eddy Simula- 
tion, small-scale turbulence is filtered out from the Navier-Stokes equations, and a 
model is used to evaluate small scales. The resulting filtered Navier-Stokes equation 
is solved for the large-scale motion, which is responsible for most of momentum and 
energy transport. The large scale of motion is highly dependent on the flow condi- 
tions and geometries under consideration. The small scales are computed from the 
turbulence model known as the subgridscale model, which in turn influence the large 
eddies. Since the small-scale eddies are more or less universal and homogeneous, it 
is postulated that the subgridscale model would be applicable to a wide range of 
flow regimes and conditions. Recall that the turbulence models for RANS are very 
much limited on the range of applications, because they attempt to model a wide 
range of scales and the random motion of eddies with no organized behavior. 

In the following sections, descriptions of filtered Navier-Stokes equations for LES 
computations and typical subgrid models are presented. 


23.2.1 Filtered Navier-Stokes Equations 


In order to explore the concept of filtering, consider a simple example of central 
difference approximation expressed as 


af| _ fag- fi- _ F(E + Az[2) - f(E — Az/2) 


Óz |; “Ar Ar (24:1) 


This expression simply approximates the ANNS at point i as the averaged values 
of dependent variable at locations i + i and i — i . Therefore, the expression (23-1) 
can be written as 


E e+ Ana- te-a 2 [xf rede] ea 


—Az/2 


which may be interpreted as a filtering operator for which any scale smaller than 
Az is filtered out. Now, a filtered quantity is defined as 


Ra) = a [^ Hae = [^57 Gta, dt (23-3) 


where, in this case, G(z,£) = 1/A, G is called a filter function, and A is the filter 
width. Typically, a filter over the entire domain is defined such that 


Fa) = f HOGE, dt (23-4) 
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It can be shown that if G is a function of (x — £), then differentiation and filtering 
operation commute [23.1]. Now, (23-4) is written as 


JG) = | HOCE- Oak (23-5) 
Some examples of one-dimensional filter functions are: 
i, 
^ if |z — ¿| < 4/2 


1. Top hat filter, G(z — £) = 
0 otherwise 


2. Gaussian filter, G(z — £) = ( )" exp [- cub 


TA? A? 
The concept of filtering is extended to three dimensions by the following 
Fa) = f FEGE- &)'& (23-6) 


where z; = x,y,z and & = 6,76 

Now, define any fluctuation at scales smaller than grid scale A as the subgrid 
scale (SGS) or unresolved quantity, and denote it by a prime. As defined previously 
by (23-4), the filtered or resolved quantity is denoted by an overbar. Thus, 


f=f+f (23-7) 


With the definitions of the resolved and unresolved quantities completed, the 
Navier-Stokes equations are now modified to yield the Filtered Navier-Stokes (FNS) 
equations. The FNS equations govern the evolution of large-scale eddies. It will be 
given for an incompressible flow initially, and, subsequently it will be extended to 
compressible flows. 

After the application of a filtering process, the incompressible Navier-Stokes 
equations become 














dn 
c0 (23-8) 

ai; 0, 1 Op Ou mn 

QU jn uc ram T a5 os) 


where 7j; is the subgridscale stress term which represents the effect of small scales. 
The subgrid scale stress is given by 


Tij = uiu; — uiu; (23-10) 


which, with the substitution of 
üi—ui-u; (23-11) 


142 Chapter 23 


can be written as 


Tj = UU; = Uü; = tuj = uj i. utu 
= Ly+Cy+ Ry (23-12) 


The terms Lij, Cij, and Rj; are defined as follows: 


Lj = üjüj —ü;iscalled Leonard stress, 
Cj = —(uuj — uju) is called cross term stress, and 
R4 = —utuj is called the subgridscale Reynolds stress 


Observe that the Leonard stress involves only the resolved quantities, and therefore 
it can be explicitly computed. Furthermore, note that this term represents the 
interaction of resolved scales which contribute and affect subgrid scales. The cross 
term stress and the SGS Reynolds stress involve unresolved quantities and must be 
modeled. The cross stress term represents the interaction of resolved and unresolved 
scales, whereas the SGS Reynolds stress represents the interaction of unresolved 
Scales. 

It is important to note that a filtered quantity represented by an overbar, that is, 
Í , is different from the averaging process used in the RANS equation, in particular 
LFF. 

It is common to rearrange Equation (23-9) and rewrite it as 

Oe Fa ia) = ve z B (p = pna) bij 
py 29 | 0 (ny - L &j) 
02;02; Oz; 4 3 EW 
1 Opt i; 0j 


^ pz; = ”Ox;02; T Oz; 








(23-13) 
where a modified pressure is defined as 


1 
p' =P- gpna (23-14) 


and i 
Tj = hj 37405 (23-15) 


pee t 
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Extension of the FNS equations to compressible flow is straightforward except 
that Favre-filtering is used. As in the case of Favre- (or mass-) averaged Navier- 
Stokes equations, this approach is taken in order to prevent introduction of sub- 
gridscale terms into the continuity equation. A Favre-filtered quantity is defined 
as 


f= e (23-16) 
Now, the Favre-filtered Navier-Stokes equations are written as 
op 


Op _ Oi | On 


dx; Oz; Oz; (23416) 


l: PEN CAESA 
a, (Pies) + Bn, Puts) + 
As in the case of incompressible flow, Equation (23-18) is rearranged as 


ð. Qi Op | ln. Of; Ory 1 OTK 
ay Pu) + dz; UU TA xi ôy 3 Oz; in dr; Ox; 36; 








or 


ais Oe pts Qo 1 OF =( 1 ) 
a 90 + Ba, Puts) i m (s = gm) ĉj = s, + as Vu 3b) (23-20) 


Define a modified pressure p* as 





Ul 
p* —B- jn (23-21) 
and 1 
T = Tij — 370 (23-22) 
Then, the momentum equation is written as 
ð =~ 0 OF Ons 
a C00 + Bz, Puts) 3 el ae bz, * Be, (23-23) 
where 
. _ 2 Óü, aù Od 
Ty = 3" 37, 9 +p (= + STA iaa Dij (23-24) 


is the filtered stress term and 


Ou; Ou; 2 du, 
Si; = Oz; += Ox; = 3 Ir, (23-25) 
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Furthermore, 
nj = p(üiü; — Wits) = pui puj/p — pt (23-26) 
is the subgrid scale stress. 

Note that at this point the subgridscale stress 7; appearing in Equation (23-20) 
is an additional unknown which must be modeled. Furthermore, the introduction of 
Tex in relations (23-21) and (23-22) and the question of how it may be computed need 
to be deliberated. These issues will be addressed shortly. For now, the goal is to 
establish the required set of filtered Navier-Stokes equations for LES of compressible 
flows. To complete the system of equations, consider the filtered energy equation 
expressed as 


ð, à m à T a 
By (Pet) + or [(pe: + pru] = B», | bo + s, 4 Siyuy) (23-27) 
where 1 1 
pe, = pë + j^ +4 0°) - 57H (23-28) 


With the assumption of perfect gas, the total energy given by relation (23-28) 
can be written as 


2 zw a E eae 1 
pe, = pT + zE T9 +H) - aT (23-29) 
which can be rearranged as 
zn fp Ek 1, 42.22 2 i: 
pà = cp ( 25) + Plt + 0° + 17^) (23-30) 


Consistent with the definition of modified pressure, define a modified tempera- 
ture as 


Papo 23-31 
2p ( ) 
Subsequently i 
pe; = cypT* + PACA T4 van (23-32) 


Consider also the equation of state for a perfect gas which, in terms of the filtered 
quantities, is written as 
p = pRT (23-33) 
In order to write the equation of state given by (23-33) in terms of modified 
pressure and modified temperature, consider the following 


. 1 enm 1 Tkk . Tkk 
Rcg ART SO cns apo E pues 
P - 37T& pR 3775 3? +5 z’ 


m(t- 25) c (R-1) 
pR (r z) + Tkk (az 3 
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or 


R 1 
+ = pRT* (zz E 3) 23-34 
p =p 4 25 3 Tkk (23-34) 


Now the energy equation given by (23-27) is written in terms of the modified 
pressure. Consider the second term given by (pe, + p)u; and rewrite it as 


(pev pui = (pec pui — (pe + p*)ü + (p& + p*)üi 








= —Qit (p& t p')ü (23-35) 
where 
Qi = —(pe + p)ui + (pe: + p*)üi (23-36) 
is the subgrid heat flux. Thus, the energy ee is now written as 
0,. On P MS e oT 
a, (Pet) + on, [(p& + p*)ū;] = um z ka +o (g Huj Sy) (23-37) 


For a relatively high Reynolds wm nion the following is introduced in 
the momentum and energy equations, respectively. 


LS = nS; (23-38) 
and . 
huj Sy = pi; Sy (23-39) 


Thus, the system of equations composed of filtered continuity, momentum, and 
energy equations is written as 


a NM 

at on, (pü;) = 0 (23-40) 

[o] prs Or ad e 

By Pts) + 5 Fa Pi + p* 3) = 7 +9 Lin Sij) (23-41) 
0; 





8 ; Dp 
a 09 + ac 2 L6 pi] = ae o,a (gr ) + (uty Sy) (23-42) 


These equations are now expanded in Cartesian coordinate and are written in a 
flux vector formulation by defining 


p pü 
pü pi? + p* 
Q=] z (23-43) E= | pao (23-44) 
pw pui 
pé (pe, + p*)ü 
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po pi 
poi pi 
F = | pù? + p* (23-45) G = | pi 
pow pi 


(23-46) 


(pe, + p*)9 (pe, -+ p*)w 
0 
Toe + Ub S 
Ey = ie T H Szy (23-47) 
$t Sox 
Q. + ore +p (aS TU Td S) 
0 
TE 4 H Ss 
F, = | Tw + Sy (23-48) 
Tye + HSyz 
Qy + KOTE + u (a8, +08, + Ü Ëp) 
y Oy T Heys yy yz 
0 
Ta + pS re 
Gy = | Tay + Ls (23-49) 
TH uS. 
Q.+ xeT* +p (üS,, t ü Sy + Ü S4) 
Now the flux vector formulation, similar to that of the Navier-Stokes equation 
given by (14-1), is written as 
8Q BE, OF 0G _ OE, 4+ 0G, 
Ot Or Oy Oz Ox Oz 
The subgridscale terms 7;; and Q; are typically expressed in terms of the eddy 


viscosity and eddy diffusivity similar to that in RANS equations and is presented 
next. 











(23-50) 


23.2.2 Subgridscale Models 


At this point, certain parallelisms between LES and RANS can be realized. 
Since a relatively strong background on RANS and turbulence models has already 
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been established, continuous reference to RANS and turbulence models will be made 
as options for LES and subgridscale models are explored. 

As in the case of RANS where turbulence quantities are written in terms of 
eddy viscosity and eddy diffusivity, the subgridscale terms are also expressed in 
terms of eddy viscosity and eddy diffusivity. For most applications in either RANS 
or LES, a model for eddy viscosity is developed, and, subsequently, eddy diffusivity 
is determined by the introduction of a turbulent Prandtl number. 

Subgridscale models similar to turbulence models vary in sophistication from 
simple algebraic models to one-equation models to multi-equation models. 

Several concepts with regard to small scales of turbulence were identified in 
Chapter One which are relevant to the development of SGS models and are re- 
viewed at this point. The small scales are more uniform and isotropic than large 
scales. Therefore, modeling of small scales by simple algebraic expression tends to 
represent the physics fairly well and is applicable to a wider range of flow condi- 
tions. Furthermore, subgridscale stress contributes a fraction of total stress, and, 
therefore, small errors in modeling do not substantially effect the overall accuracy. 

In addition to assumptions of isotropic and uniformity of small scales, the as- 
sumption of equilibrium is also imposed for most applications. The assumption of 
equilibrium is referred to as a condition where the energy contained in large ed- 
dies is transferred to smaller eddies, and, ultimately, all are dissipated into heat at 
the level of smallest eddies (subgridscale) due to molecular viscosity. Although, on 
the average, the transfer of energy takes place from resolved scales to unresolved 
scales, the reverse can occur (known as backscatter), which is not accounted for 
when equilibrium assumption is used. 

Based on these discussions, algebraic models appear to be good candidates for 
subgridscale terms. Besides their simplicity, they require a minimum amount of 
computational effort. In the following sections, algebraic SGS models are intro- 
duced. 


23.2.2.1 Eddy Viscosity 


The subgridscale stress is related to the large-scale strain rate tensor by introduction 
of eddy viscosity and is written in a similar form as the laminar stress. For an 
incompressible flow, the subgridscale stress term is written as 


1 z 
T = Tij = 3 Tkk bij = 2 Sij (23-51) 


where 





S; = ; (= + S) (23-52) 
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The subgridscale stress term for compressible flow is written as 
Ty = Tij — ira ó = Du S (23-53) 
and Si is given by : . 
Sy = Be: + os x (23-54) 


The subgrid heat flux is written in terms of eddy diffusivity and, subsequently, 
in terms of turbulent Prandtl number as follows. 





(23-55) 


Observe that, just as in the RANS equations, the eddy conductivity k, is replaced 
by terms of eddy viscosity and turbulent Prandtl number. Again, just as in the 
RANS equations, several models are available to determine v. 


23.2.2.1.1 Smagorinsky Model 


The first subgridscale stress model was developed by Smagorinsky [23.2] in the 
1960's. It is a simple model which provides reasonable results for some applications. 
This model is still commonly used in LES due to its simplicity. Several modifications 
have been proposed to improve its prediction capabilities. One such attempt is the 
dynamic model which will be discussed in the next section. 

The Smagorinsky model is developed based on the assumption of equilibrium. 
The eddy viscosity is written to be proportional to a length scale £. 

The subgridscale stress is written in terms of eddy viscosity as 





1 z 
Tij — 3 Tkkôij = 25, Sij (23-56) 
where, for an incompressible flow, the large-scale strain rate is given by 
= lfí0ü , Ou; 
m — 23-57 
Sij 2 (s m 24) ( ) 
and in terms of the velocity scale qscs as follows 
Wc lasas (23-58) 
The length scale is taken as the filter width A, and the velocity scale gses ~ £ (Sl, 
where 
|5| = 254 Sy (23-59) 
Finally, 


v ~ A’ S] (23-60) 
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which is written as an equation by introduction of the Smagorinsky constant as 
— 2 AG 
n = c A*|S| (23-61) 


The Smagorinsky constant can be approximated in terms of Kolmogorov con- 
stant cy by the following 
»]j3 VM 
C, 2 — (5 cx) 
T 


2 
where, if c, = 1.4, then c, = 0.18. This value of Smagorinsky constant is not 
universal. In fact, similar to the constants used in the algebraic turbulence models, 
it’s value would vary. However, for most applications, a constant in the range of 
0.1 « c, « 0.24 have been used. 
The Smagorinsky model and its application are simple and provide reasonable 
results for some flows. However, there are several drawbacks to the model. 


e First, as mentioned above, the value of Smagorinsky constant is not universal, 
and its value is different for different types of flows. 


e Second, due to the imposed equilibrium assumption, the subgridscale stress 
always decreases the energy of flow, and, therefore, backscatter which may be 
present in the flow is omitted. 


e Third, it may be necessary to reduce the model constant in the near wall 
region. 


e Fourth, since, for most applications, grid clustering will be used, specification 
of filter width A may be difficult. Among several options available, one may 
consider A = (Az? + Ay? + Az?) ? or A = (AzAyAz)!5. 


An improvement to the model can be made where the value of constant c, would 
vary locally, thereby reducing or increasing eddy viscosity in locations as required. 
This concept is the basis of the dynamic model which is provided next. 


23.2.2.1.2 Dynamic Model 


Several potential problems with the original Smagorinsky model may develop in 
LES of near wall flows, shear flows, or in transition regions as discussed in the 
previous section. These difficulties are due to the selection of a constant value for 
Smagorinsky coefficient. From physical observations, it is known that, for example, 
backscatter occurs in some flows which in fact can be significant. Furthermore, 
the subgridscale stress must possess assymptotic behavior near the wall and vanish 
at the wall To account for these physical considerations and to overcome the 
difficulties associated with a constant value of c,, a dynamic model is proposed 
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by Germano [23.3]. In this model, the coefficient c, is allowed to vary locally to 
adjust the level of eddy viscosity in the flowfield. The procedure is initiated by 
specifying a value of c,, and subsequently after each time level or in practice after 
several time levels, the flowfield is examined, and a new value of c, at each point 
is determined. Therefore, the value of the coefficient c, is continuously modified 
as solution proceeds in time. The scheme is similar in principle to adaptive grid, 
where the grid system evolves as the solution proceeds in time. 

The dynamic model is developed by the introduction of a test filter with a length 
scale larger than that of the length scale of the original filter, typically by a factor 
of two. Stresses in this range are referred to as subtest-scale (STS) stresses and are 
modeled similar to the SGS stresses. The test filter is defined by G, and a tophat 
*^" is used to identify the associated test-filtered quantities. After the test filter is 
applied to the FNS equations, the STS stress is developed and is given by 


Fy = tt; — uiu; (23-62) 
and the resolved stress is 
Ly = ty; — a; (23-63) 


Recall that the SGS stress as given by (23-10) is 

nj = üjü; — HG (23-64) 
When the test filter is applied to (23-64), one obtains 

fy = dii — uj (23-65) 


It is then obvious that one can write 


a 


Li iun Tig = =f (23-66) 


The resolved stress Li can be computed explicitly from the resolved velocity 
field, whereas the SGS stress and STS stress require modeling. 

In the following, the dynamic model is used in conjunction with the Smagorinsky 
model. In fact, the procedure can be applied to any eddy viscosity model. Now, 
the Smagorinsky model applied to SGS and STS stresses is written as 


gr olx = 
Tij — 3 Tkk OFF = 26 ij (23-67) 


and " 
Tg — 3 Ty. bij = 2c D (23-68) 
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where 
Qi = A? I| S; (23-69) 
By = A? ls | S; (23-70) 


Subtracting (23-67) from (23-68), one obtains 


"E ee _ 
Tg — Ty — 3 (7. 1 fu) ó; = 2c B; — ai; (23-71) 
or i 
Lg — 3 Lin by = 2c By — 2603 (23-72) 
Typically, c is removed from the filtering operation such that (23-72) is written 
as : P 
Lij — 3 Lin ij = 2c (By — à) = 2c Mi; (23-73) 
where 


My = By — à (23-74) 
Equation (23-73) yields five independent equations with only one unknown c, 


that is, the system is overdetermined. One approach to overcome this problem is 
proposed by Germano, in which Equation (23-73) is contracted by Sj, resulting in 





Êy Sy = 2c My S (23-75) 
from which is 
— 1 Ly Si 

cms My Ss (23-76) 


Another approach is proposed by Lilly [23.4], in which the sum of the squares 
of the residual is minimized. The result is 





(23-77) 


This formulation is particularly attractive because the formulation given by (23-76) 
may cause a difficulty in that the denominator could become locally zero or very 
small such that numerical instability within the solution develops. This difficulty 
in formulation (23-76) can be removed by the assuniption that c is a function of y 
and t only. This assumption allows the introduction of an averaging scheme taken 
over the z and z directions, that is, parallel to the wall. This averaging scheme 
(indicated by < >) also addresses another difficulty which may be encountered 
when either (23-76) or (23-77) is used. It has been observed that the value of c 
could vary significantly over the domain and may possess negative values over large 
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regions. This may also become a source of instability in the solution. Though a 
certain value of negative c is desireable in the solution, large values could create 
serious problems. The negative values of c are interpreted as allowing backscatter 
in the solution. Introduction of the averaging procedure discussed above tends to 
eliminate the problem of excessive negative values of c. However, note that, with 
the introduction of the averaging process over the domain, the model can no longer 
be considered strictly local. 

Now the formulations of c are written in terms of the averaging scheme as follows. 


_l< Ly Si > 
ena MEUS (23-78) 
and y 
BE Rae (23-79) 


rur Mè > 

In closing this section, it is concluded that the dynamic model overcomes some 

of the problems associated with eddy viscosity models with constant coefficient. In 

the process, however, the dynamic model also creates some difficulties as identified 

above. However, some simple procedures were introduced by which these difficulties 
can be removed. 


23.3 Direct Numerical Simulation 


The third category of solution scheme for laminar, transitional, and turbulent 
flows is Direct Numerical Simulation. This scheme essentially solves the time depen- 
dent Navier-Stokes equations and attempts to accommodate all scales of turbulence. 
There are two major grand challenges that need to be resolved in order for DNS 
to be a routine computational scheme for design and analysis purposes. The first 
major obstacle is the limitation of the available resources. Present day computers 
do not have the memory, operation speed, and data transfer capabilities to conduct 
a DNS of complex configuration at high Reynolds number. 

It was shown in Chapter 20 that the number of grid points required for DNS 
is proportional to Re?/*, For example, a DNS of a complete aircraft will require 
about one extaflop, that is 10!? flops. The projected computer performance is 
approximately 10° (operations per second) to be achieved in about 2010. One 
approach to gain higher performance and increase memory is to perform DNS on 
multi-processor, parallel computers. In this approach, the domain of solution is 
divided among several processors, where information is transferred between the 
processors. And, of course, consideration of data transfer and the transfer speed is 
crucial. However, multi-processor parallel computations have shown great potential 
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in performance of DNS for turbulent flows over large domain at high Reynolds 
numbers. 

The second challenge is the development of a solution algorithm which is free 
of significant numerical error. This aspect of the problem includes the grid system, 
numerical scheme, and boundary conditions. 

The grid system must be developed by a higher-order grid generation scheme 
such as to provide a high quality grid, that is, a grid for which there is continuity in 
the curvature and the metrics obtained by high-order approximations, for example, 
fourth-order or higher. An example of such grid generation is the fourth-order 
elliptic grid generation scheme. 

The development of numerical schemes must include higher-order approximation 
of the derivatives. The fourth-order Runge-Kutta scheme can be used for time 
integration, whereas typically sixth-order compact finite difference approximation 
is used for the convection and diffusion terms. The third element of the solution 
algorithm which is crucial is the treatment of boundary conditions. The sponge layer 
type boundary treatment is usually implemented for the inflow/outflow boundaries, 
and the standard no-slip boundary condition is applied at the solid surface. 

Some examples of DNS of transitional/turbulent flows which have been reported 
in the literature are provided below. As the technology matures, the envelope and 
the range of applications will increase to include higher Reynolds numbers and more 
complex geometries. The following list is a very limited example of DNS and by no 
means is complete. 


* Compressible free shear flow. Several computations are performed at convec- 
tive Mach number (M.) of 0.38 (M, = 1.5, M; = 1.5), M. of 0.4 (M; = 2.0, 
M2 = 1.2), and M, of 0.8 (Mi = 4.0, M; = 2.4). The Reynolds number range 
was between 100 and 500, where Re — pi(ui — u3)ó6,/ui1. The subscripts 1 and 
2 refer to high and low speed streams, and 6, is the initial vorticity thickness 
[23.5]. 


* Transitional flow at Mach 0.5 and He; of 160,000 over a Joukowsky airfoil. 
The effect of suction/blowing is considered as à means for controlled transition 
[23.6]. 


* Flow transition at Mach 1.6 over a Joukowsky airfoil The Reynolds number 
based on the half-thickness of the airfoil A was Re, = 15000 [23.7]. 


* Compressible flat plate boundary layer flow at Mach 2.25. The Reynolds 
number based on the inlet conditions was 635000/in [23.8]. 


e Fully turbulent channel flow. Effect of wall injection and suction was investi- 
gated. The Reynolds number based on the channel half width and averaged 
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friction velocity was 150 [23.9]. 


e Supersonic turbulent boundary-layer flow over a flat plate and a compression 
ramp. The compressible, turbulent flowfield over a flat plate at Mach numbers 
of 3, 4.5, and 6, and at momentum thickness Reynolds numbers of 3015, 
2618, and 2652 were investigated. In addition, flowfield over an 18 degree 
compression ramp at Mach 3 and Reynolds number of 5000 was considered 
[23.11]. 


pe 


APPENDIX J: 


Transformation of Turbulence Models from 


Physical Space to Computational Space 





The details of the coordinate transformation for turbulence models are provided 
in this appendix. Since the terms for most models in each category, e.g., one- 
equation models, are similar, mathematical steps will be carried out for a typical 
model only. Extension to similar turbulence models is straightforward. 


J.1 Baldwin-Barth Turbulence Model 


To better manage the mathematical work, each term in the equation will be 
considered separately. Again, using the notation R for vRer, the left-hand side of 
(21-86) is dR/dt, which is expanded in Cartesian coordinates to yield 


dR ðR y OR ƏR OR 
da ber VVR-Uy tug tis 
Now, using the Cartesian derivatives (21-88) and (21-89), the convective term be- 
comes 
dR | OR _ pons OR my LL OR 
Ur Dy T dE ^ Gy) TY CA: TOR 
OR OR OR ve 
= (&zu + &,v) OE + (neu + Ty) ôn -UX Y 


The diffusion term in Equation (21-86) involves the Laplacian V?R and Vu - VR. 
Each term is treated separately, and details are as follow. 
Using the expressions given by (21-90) and (21-91), one has 


8R ms OR "OR WR 


?R m 
V ag ag ME t Eg + 


= Oar 
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BE OE | BE On 


T 06, OR | On, OR 06. OR | On. OR 
On 06 n On 


eR OR  ,0R 
LOL E + yw Eon ur 


& 96, OR , Ony OR d& OR , On, OR 
8€ OF | OE On "10g O€ 0g An 


Rearrange terms, 


OR 


OR R 
VR = (+g) de + Gere + im) gens + (nf + n) att 


E df: OR , Ons OR , ( OR q 9n. OR 
"\ GE OE © 0€ On On eet ‘On On 


+6, (25 OR , 9n, ORY |, (86, OR | 0n, OR 


tl 


(€ t eee m + 2(£m. + Ene e ran + (n2 + 2s an i 


Ee, x BE, A) aR 
+(65 "get 8g ac aE t 8€ 

z z R 

+ (eat + nee + Ty 3) » 


Using previously defined expressions given by (11-210d), (11-211d), and (11- 
212e) repeated here for convenience, and additional definitions, the Laplacian be- 
comes 


OR OR OR OR 


u ZR 
2 
V"R- Oa Ber + 326569, t ae 9 GE t+ 9 Gy 
where 
a = E+E 
b, = m+n 
cs = Ene + Eyty 
= oe uS 4 28 4 
n= Gg tm a, + ge thg 
= gz p 2n, g Ith on 
g = ex + "n + éve Ön 
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The second term Vr - VR is first expanded as 


àv, OR Ou, OR 


Vw YR = ĈA IR ds OR 
"Nom Oy Oy 


Now, using the Cartesian derivatives, one has 


oy, oy, OR OR ov, an, OR | OR 
(ext + ng) (ex 8€ ate Ne Sp 3E +(e% ðE Th 2a) G 8€ + y =) = 


Eon DR OR re ov, OR | nat OY OR | 20% OR 
"8p of 0E On | H On OF ME dn On 
Ou, OR ov, OR , Ov, OR | ôn OR 
9v V VeL uiv t n t 
ur d thag an +O BE Gy On 
Ov, OR Ov, ôR 
_ 2 : evt evt 


ov, OR ðv: ôR 
F (Ens + &yny) ôn 0€ + (nz + n) ôn On 


Lan OR , 8w BR _ Ow DRY , , Ou OR 
~ OE OE OE ðn On OE) “On ðn 


Now, consider the Baldwin-Barth turbulence model given by (21-87). The trans- 
formation of the convective term is exactly the same as the previous case, that is 


= OR OR OR OR 
Patena gh pe og 


The Laplacian in the diffusion term is first expanded in Cartesian coordinates 
and expressed as 


OR OR ð (dR 8 (OR 0A - BB 
2 OH GH O0f[On|, O[0H 
VRTM +e dx (=) (#)- Or By dy 
where 
OR OR 


Now relations (21-88) and (21-89) are used to provide 

ðA ðA ðB ôB 
2 L3] — ———— eee —— 
VR = EE MET + rae MET 


ð 0A ðB 
= E73 + &5¢ + EDI + On 
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The second expression in the diffusion term V - V, VR is written as 


V.uVR = V- CI i) 


Oz Oy 
_ ô (OR), 8 (ƏR) _ ƏC , aD 
~ Or V az Oy V Oy] ax" dy 
where 
OR OR 
C= ae and D= "ay 
In the generalized n one i 
ðc ôD 
V-nVR = e TU» T 
E 9 (eee OR ic d cbe on OR 
= BE t | Sz 8E Nra an "On t | Sz 8E Nz 8n 
" m on Paci OR " LA UM OR 
OE t| Sy ET WA an "ay t | Sy 8€ Wa an 
m op 0p 
where 


OR 
Q&Q = h (e x + Ne x) 
B Lu G + Wat 2) 


Finally, the following definitions are introduced which will be used in the formulation 
of FDE. 


Il 


lu pian € Vna 
dim Gla T 


"TES 88 oa, , 2B 


Rew 9. 


wager A) 8B ðA OB 
MUT VA eae t frag + Gy + Gy 


The last term of the Baldwin-Barth turbulence model includes a production 


term given by 
(ce fo — ca)y Ca Dı Di vRer S 
where S in the computational space is given by (21-98). 





APPENDIX K: 


The Transport Equation for the 
Turbulence Kinetic Energy 





It is well established that turbulence in a flowfield is generated due to a variety 
of factors such as, for example, surface roughness, and it is convected and dissi- 
pated throughout the domain. To represent the above-mentioned physical aspect 
of turbulence in the development of a turbulence model, a transport equation is 
developed in terms of the turbulence kinetic energy. This transport equation which 
is derived from the Navier-Stokes equation is referred to as the turbulence kinetic 
energy equation (TKE). The derivation of turbulence kinetic energy and the inter- 
pretation of terms in the equation are the focus of this section. The derivation of 
the equation will be performed in tensor notation. Recall that a repeated index in 
tensor notation indicates summation over that index. 

The contintuity equation and the i-component of the momentum equation ex- 
pressed in tensor notation are 


Op. Oa 
at + Óz; (puj) =0 (K-1) 
"e 0 ə dp ð 
i yu, OM occ OR, Ong E 
p ( ot t uj 22) Oz; T Oz; (K 2) 
where à 8 ra 
ep eee) a6. 4 Ok E 
a (5 t 24) $534 A os) 


The Reynolds averaged continuity and the i-component of momentum equations 
are given by 


SE REN RE 
z T Ba, 0) + Ba, V) =0 (K-4) 
and 
_ 04; ðu; PS T Ou; 
Pup ug, T VERUS 


160 Appendix K 








=- R + "i (K-5) 


To proceed with the derivation of a TKE equation, multiply the continuity 
equation by 1/2 u? to obtain 








15,|0p ô 
2" lara ow | =0 (K-6) 
Multiply the i-component of the momentum equation by u; to obtain 
Ou; Ou; Op Ónj 
; ; oy d up K- 
pu (5 was) "ass ir, den 


Note that the second term in Equation (K-7) can be rearranged as 


Vui DC os LI EDI 
Peis On, TP A2 Oz; 


Thus, Equation (K-7) is written as 


uj on 
pui gi 3 Pi am i (u?) = -u =— = P tu ot Ox; (K-8) 
Add Equation (K-8) to Equation (K-1) to obtain 
ð u? 8 2\ _ Op , On; 
at (d Dex (ou) = -u 5 + "i Oz; (K-9) 


Now substitute the sum of average and fluctuating quantities for the instanta- 
neous quantities, so that Equation (K-9) becomes 


ð |} : - ð [1,. » 7 n2 
5 5 OHA) Gcr ul? | n [5 0-6 Gy + uh) (+ ut) 
Ot |2 Oz; 2 j 7 
- (tat ot) a (p p) + (i ui) 5— Uv Tt 5) 
Taking the time average following the rules set " (17-2), one obtains 
18 = ——À — 
za; (Bt put gut 2u pul) 


l 89 p.56. _—F aes 
+ 2 Ba, [P (ast? + au?) + uspu, + 2u;üjp'ui 
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+ puju + 2pü;uju; + üj pu, + pula? + 2u, ufu] 


zi Op , OP = OT „OTi; 
Sg ee: Peek o, K-10 
“On, us Ti rs T ar ( ) 


Now multiply the momentum equation given by Equation (K-5) by ŭ;. 


.([.9ü Ou fn, ao OU 
a (a FE ere) en (pa, + 7 bz, 


+ izz (z uuj) 4 M (ipa) + T ( ufu) 
- -us + ag (K-11) 


ion. IO pam. 18 y ii 
aUa taa OM) tx (Put) +a Ge 
1 8 (Ot; 


1 ey SN ew SI 
2 da, [Pajuk + aul + pull + puut] 


Ax; ‘Ər (K-12) 
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The fifth and sixth terms can be combined as follows 


1 ð __. 0&6 1,0 ,__ 
2 Oz; (puja ) OP = 3' 37; (pu;) 
The remaining terms can be combined in a similar fashion, and Equation (K-12) 
written as 
1_,0p 108 
iU tH OM) +55 OH) rr 
T Ou 1_, 0 
+ i 


9 ig (ju) + Pung 5 + 3a * Oz; 


; (4) 


1 LLL. 
20 [pau + a; pul pil + pujut + püjut | 
Tj 


Op Or; 
= -v— uae - 
s uri + “Oz; (K 13) 


Observe that the sum of the first and fifth terms is zero, due to continuity, and one 
can write 


Op’ 8 Emi 
—-u— = —— (pul x -14 
ul bz, Ce (rui) +p dn, (K-14) 
Now define the turbulent kinetic energy k as 
1 
k=5 (uf) (K-15) 
Finally, Equation (K-13) is written as 
K+ a (Pueri) = -a a Gu ul) 
- Ot; 
— = (puj k) "AL = (ug) - puu Ua oz; 
- gui vu 2s b pq P ey euo 198 (K-16) 
53 Ax "5 os; 


This is the turbulence kinetic energy equation. 
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Recall that the shear stress is given by 








= Ou; Ou; 2 Ou, 
Ce (ss 5 Ox; by 3" Oz, qe 
Furthermore, one can write 
Ori, à Ovi 
pS = —— (ulrl) pt Sot 4 
VAL = Bay (ulti) - mh s, (K-18) 


Now, an alternate form of the TKE equation can be written by substitution of 
(K-17) and (K-18) into (K-16), as follows 
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The interpretation of the terms in the TKE equation given by (K-19) is as 
follows. 


I. The rate of change of TKE which is composed of a local change and a con- 
vective term is 


x (pk) aa — (pu Ñ; k + p'u} gj k) (K-20) 
II. The diffusion term includes the following 


a.u ul) — 3 a; OR k) — ie (a; pk) (K-21) 


and represents diffusion of turbulence due to molecular transport processes. 
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III The production term is 


——7 Ou; Ou; 
— puiu; ira - (y Pg pu; + puiu pu) a (K-22) 


and it represents the rate of transfer of kinetic energy of the mean flow to 
turbulence. 


IV. The work done by the turbulent viscous stresses is given by 








ie fu g (Fe + 54) - bij : E } (K-23) 
V. The turbulence viscous dissipation is 
ERE w 
and finally 
VI. The pressure dilatation is 
yet (K-25) 


Typically the fluctuating correlations involving density fluctuations are dropped. 
Furthermore, the Boussinesq assumption is invoked, and the terms defined by (K- 
20) through (K-25) are written as 


I. Rate of change of TKE is 
ô ,. à uu. 
a; (PR) + 3a; (pu; k) (K-26) 
Note that this expression can be modified as follows 


_ Ok Op _ Ok Ó ,.- 
2 GE + p NET bt MU dr (pu;) 


dk 
=p (G+ je) =F (K-27) 


where continuity is used to set the bracket term equal to zero. 
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IL The diffusion term is written as 


ô [fq Ok 
Oz; (4 =) (K-28) 
III. The production is denoted by P, and is given by 
à Oü; 
Pr = —pujul —— = Tiy — (K-29) 
v Oz; z Oz; 
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or 
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IV. The work done by the turbulent viscous stress is 


à Ok 
az (^35) iss 
V. The turbulent viscous dissipation, which is typically referred to as dissipation, 


is 
Ou; Ou; = 
—p 32,02, a —pe (K-33) 


VI. The pressure dilatation term is usually dropped. Finally, the TKE equation 
is expressed as follows 
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